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Abstract 

We deal with the problem of bridging the gap between two scales in neuronal 
modeling. At the first (microscopic) scale, neurons are considered individually 
and their behavior described by stochastic differential equations that govern 
the time variations of their membrane potentials. They are coupled by synap- 
tic connections acting on their resulting activity, a nonlinear function of their 
membrane potential. At the second (mesoscopic) scale, interacting populations 
of neurons are described individually by similar equations. The equations de- 
scribing the dynamical and the stationary mean field behaviors are considered 
as functional equations on a set of stochastic processes. Using this new point of 
view allows us to prove that these equations are well-posed on any finite time 
interval and to provide a constructive method for effectively computing their 
unique solution. This method is proved to converge to the unique solution and 
we characterize its complexity and convergence rate. We also provide partial 
results for the stationary problem on infinite time intervals. These results shed 
some new light on such neural mass models as the one of Jansen and Rit j2?| : 
their dynamics appears as a coarse approximation of the much richer dynamics 
that emerges from our analysis. Our numerical experiments confirm that the 
framework we propose and the numerical methods we derive from it provide a 
new and powerful tool for the exploration of neural behaviors at different scales. 
Keywords: mean-field analysis, stochastic processes, stochastic differential 
equations, stochastic networks, stochastic functional equations, random con- 
nectivities, multi populations networks, neural mass models. 
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1 Introduction 



Modeling neural activity at scales integrating the effect of thousands of neurons 
is of central importance for several reasons. First, most imaging techniques are 
not able to measure individual neuron activity ("microscopic" scale), but are 
instead measuring mesoscopic effects resulting from the activity of several hun- 
dreds to several hundreds of thousands of neurons. Second, anatomical data 
recorded in the cortex reveal the existence of structures, such as the cortical 
columns, with a diameter of about 50//m to 1mm, containing of the order of 
one hundred to one hundred thousand neurons belonging to a few different 
species. These columns have specific functions. For example, in the visual cor- 
tex VI, they respond to preferential orientations of bar-shaped visual stimuli. In 
this case, information processing does not occur at the scale of individual neu- 
rons but rather corresponds to an activity integrating the collective dynamics of 
many interacting neurons and resulting in a mesoscopic signal. The description 
of this collective dynamics requires models which are different from individual 
neurons models. In particular, if the accurate description of one neuron requires 
"m" parameters (such as sodium, potassium, calcium conductances, membrane 
capacitance, etc.), it is not necessarily true that an accurate mesoscopic de- 
scription of an assembly of N neurons requires Nm parameters. Indeed, when 
N is large enough averaging effects appear, and the collective dynamics is well 
described by an effective mean-field, summarizing the effect of the interactions 
of a neuron with the other neurons, and depending on a few effective control pa- 
rameters. This vision, inherited from statistical physics requires that the space 
scale be large enough to include a large number of microscopic components (here 
neurons) and small enough so that the region considered is homogeneous. This 
is in effect for instance the case of cortical columns. 

However, obtaining the evolution equations of the effective mean-field from 
microscopic dynamics is far from being evident . In simple physical models this 
can be achieved via the law of large numbers and the central limit theorem, 
provided that time correlations decrease sufficiently fast. This type of approach 
has been generalized to such fields as quantum field theory or non equilibrium 
statistical mechanics. To the best of our knowledge, the idea of applying mean- 
field methods to neural networks dates back to Amari [1, Q . In his approach, 
the author uses an assumption that he called the "local chaos hypothesis", 
reminiscent of Boltzmann's "molecular chaos hypothesis" , that postulates the 
vanishing of individual correlations between neurons, when the number TV of 
neurons tends to infinity. Later on, Crisanti, Sompolinsky and coworkers 34 1 
used a dynamic mean-field approach to conjecture the existence of chaos in an 
homogeneous neural network with random independent synaptic weights. This 
approach was formerly developed by Sompolinsky and coworkers for spin-glasses 
J3, 17, 3H, where complex effects such as aging or coexistence of a diverging 



number of metastable states, renders the mean-field analysis delicate in the long 
time limit 

On the opposite, these effects do not appear in the neural network considered 
in 34 1 because the synaptic weights are independent [l(j (and especially non 
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symmetric, in opposition to spin glasses). In this case, the Amari approach 
and the dynamic mean-field approach lead to the same mean-field equations. 
Later on, the mean-field equations derived by Sompolinsky and Zippelius 13511 
for spin-glasses were rigorously obtained by Ben Arous and Guionnet 0, [H, l24j . 
The application of their method to a discrete time version of the neural network 
considered in 34 1 and in [3l[ was done by Moynot and Samuelides [32[ . 

Mean-field methods are often used in the neural network community but 
there are only a few rigorous results using the dynamic mean-field method. The 
main advantage of dynamic mean-field techniques is that they allow one to con- 
sider neural networks where synaptic weights are random (and independent). 
The mean-field approach allows one to state general and generic results about 
the dynamics as a function of the statistical parameters controlling the prob- 
ability distribution of the synaptic weights [33j . It does not only provide the 
evolution of the mean activity of the network but, because it is an equation 
on the law of the mean-field, it also provides informations on the fluctuations 
around the mean and their correlations. These correlations are of crucial impor- 
tance as revealed in the paper by Sompolinsky and coworkers j34[. Indeed, in 
their work, the analysis of correlations allows them to discriminate between two 
distinct regimes: a dynamics with a stable fixed point and a chaotic dynamics, 
while the mean is identically zero in the two regimes. 

However, this approach has also several drawbacks explaining why it is so sel- 
dom used. First, this method uses a generating function approach that requires 
heavy computations and some "art" for obtaining the mean-field equations. Sec- 
ond, it is hard to generalize to models including several populations. Finally, 
dynamic mean-field equations are usually supposed to characterize in fine a 
stationary process. It is then natural to search for stationary solutions. This 
considerably simplifies the dynamic mean-field equations by reducing them to a 
set of differential equations (see section [5]) but the price to pay is the unavoid- 
able occurrence in the equations of a non free parameter, the initial condition, 
that can only be characterized through the investigation of the non stationary 
case. 

Hence it is not clear whether such a stationary solution exists, and, if it 
is the case, how to characterize it. To the best of our knowledge, this difficult 
question has only been investigated for neural networks in one paper by Crisanti 
and coworkers [15[ . 

Different alternative approaches have been used to get a mean-field descrip- 
tion of a given neural network and to find its solutions. In the neuroscience 
community, a static mean-field study of multi population network activity was 
developed by Treves in [36| . This author did not consider external inputs but 
incorporated dynamical synaptic currents and adaptation effects. His analy- 
sis was completed in where the authors considered a unique population of 
nonlinear oscillators subject to a noisy input current. They proved, using a 
stationary Fokker-Planck formalism, the stability of an asynchronous state in 
the network. Later on, Gerstner in (2lj built a new approach to characterize 
the mean-field dynamics for the Spike Response Model, via the introduction 
of suitable kernels propagating the collective activity of a neural population in 
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time. 

Brunei and Hakim considered a network composed of integrate-and-fire neu- 
rons connected with constant synaptic weights @. In the case of sparse con- 
nectivity, stationarity, and considering a regime where individual neurons emit 
spikes at low rate, they were able to study analytically the dynamics of the 
network and to show that the network exhibited a sharp transition between 
a stationary regime and a regime of fast collective oscillations weakly synchro- 
nized. Their approach was based on a perturbative analysis of the Fokker-Planck 
equation. A similar formalism was used in [30| which, when complemented with 
self-consistency equations, resulted in the dynamical description of the mean- 
field equations of the network, and was extended to a multi population network. 

Finally, Chizhov and Graham (l3j have recently proposed a new method 
based on a population density approach allowing to characterize the mesoscopic 
behaviour of neuron populations in conductance-based models. We shortly dis- 
cuss their approach and compare it to ours in the discussion section [6] 

In the present paper, we investigate the problem of deriving the equations of 
evolution of neural masses at mesoscopic scales from neurons dynamics, using 
a new and rigorous approach based on stochastic analysis. 

The article is organized as follows. In section [5] we derive from first princi- 
ples the equations relating the membrane potential of each of a set of neurons as 
function of the external injected current and noise and of the shapes and inten- 
sities of the postsynaptic potentials in the case where these shapes depend only 
on the post-synaptic neuron (the so-called voltage-based model) . Assuming that 
the shapes of the postsynaptic potentials can be described by linear (possibly 
time-dependent) differential equations we express the dynamics of the neurons 
as a set of stochastic differential equations. Assuming that the synaptic connec- 
tivities between neurons satisfy statistical relationship only depending on the 
population they belong to, we obtain the mean-field equations summarizing the 
interactions of the P populations in the limit where the number of neurons tend 
to infinity. These equations can be derived in several ways, either heuristically 
following the lines of Amari 0,01, Sompolinsky [H, H3, and Cessac [loll33j| . or 
rigorously as in the work of Benarous and Guionnet [J, 0| . The purpose of 
this article is not the derivation of these mean-field equations but to prove that 
they are well-posed and to provide an algorithm for computing their solution. 
Before we do this we provide the reader with two important examples of such 
mean-field equations. The first example is what we call the simple model, a 
straightforward generalization of the case studied by Amari and Sompolinsky. 
The second example is a neuronal assembly model, or neural mass model, as 
introduced by Freeman [20| and examplified in Jansen and Rit's cortical column 
model [27j] . 

In section [3] we consider the problem of solutions over a finite time interval 
[io,T]. We prove, under some mild assumptions, the existence and uniqueness 
of a solution of the dynamic mean-field equations given an initial condition at 
time to. The proof consists in showing that a nonlinear equation defined on the 
set of multidimensional Gaussian random processes defined on [to, T] has a fixed 



4 



point. We extend this proof in section 0] to the case of stationary solutions over 
the time interval [—00, T] for the simple model. Both proofs are constructive and 
provide an algorithm for computing numerically the solutions of the mean-field 
equations. 

We then study in section \5\ the complexity and the convergence rate of this 
algorithm and put it to good use: We first compare our numerical results to 
the theoretical results of Sompolinsky and coworkers IB], 34 1 . We then provide 



an example of numerical experiments in the case of two populations of neurons 
where the role of the mean-field fluctuations is emphasized. 

Along the paper we introduce several constants. To help the reader we have 
collected in table[T]the most important ones and the place where they are defined 
in the text. 



2 Mean-field equations for multi-populations neu- 
ral network models 

In this section we introduce the classical neural mass models and compute the 
related mean-field equations they satisfy in the limit of an infinite number of 
neurons. 

2.1 The general model 
2.1.1 General framework 

We consider a network composed of N neurons indexed by i € {1, . . . , N} 
belonging to P populations indexed by a £ {1, . . . , P} such as those shown 
in figure [T] Let N a be the number of neurons in population a. We have 
N = Yla=i -^a- We define the population which the neuron i, i = 1, • • • , N 
belongs to. 

Definition 2.1. The function p : {1, . . . , N} — > {1, . . . , P} associates to each 
neuron i £ {1, • • • , N}, the population a = p(i) £ {1, • • • , P}, it belongs to. 

We consider that each neuron i is described by its membrane potential Vi (t) , 
and the related instantaneous firing rate is deduced from it through a relation 
of the form Vi{t) = Si(Vi(t)) [H], [22|, where 5, is a sigmoidal function. 
A single action potential from neuron j generates a post-synaptic potential 
PSPij(u) on the postsynaptic neuron i, where u is the time elapsed after the 
spike is received. We neglect the delays due to the distance travelled down the 
axon by the spikes. 

Assuming that the post-synaptic potentials sum linearly, the average mem- 
brane potential of neuron i is 

Vi(t)= PSPij(t-t k ) + Vi(to), 

j,k,t k >t 
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Figure 1: General network considered: N neurons belonging to P populations 
are interconnected with random synaptic weights whose probability distribu- 
tions only depend upon the population indexes, see text. 

where the sum is taken over the arrival times of the spikes produced by the 
neurons j after some reference time to. The number of spikes arriving between 
t and t + dt is Vj(t)dt. Therefore we have 

Vi( t ) = J2 I PSPi J (t-s)is J ( s )ds+V l (t ) = Y, I PSP ij {t-s)S j {Vj{s))ds+V i {t ), 

„• Jto j Jto 

(1) 
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or, equivalcntly 



Vi(t) = S 1 [y^J^ PSP l3 {t - s)v 3 {s) ds + Vi{to) \ . (2) 

The PSPijS can depend on several variables in order to account for instance for 
adaptation or learning. 

We now make the simplifying assumption that the shape of the postsynap- 
tic potential PSPij only depends on the postsynaptic population, which corre- 
sponds to the voltage based models in Ermentrout's classification 



The voltage-based model 

The assumption, made in (25l |. is that the post-synaptic potential has the same 
shape no matter which presynaptic population caused it, the sign and amplitude 
may vary though. This leads to the relation 

PSP ij {t) = J ij9i {i). 

gt represents the unweighted shape (called a g-shape) of the postsynaptic po- 
tentials and Jij is the strength of the postsynaptic potentials elicited by neuron 
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j on neuron i. At this stage of the discussion, these weights are supposed to be 
deterministic. This is reflected in the notation which indicates an average 
value. From equation JT]) we have 



^J y -zA,(s) ds + Vi(t ). 



So far we have only considered the synaptic inputs to the neurons. We enrich our 
model by assuming that the neuron i receives also an external current density 
composed of a deterministic part, noted /*(£), and a stochastic part, noted rii(t), 
so that 




Vi(t) = jf gi (t - s) far JijVjis) + Ii(s) + ni {s)J ds + Vi(t ). (3) 

We assume, and this is essential for deriving the mean-field equations below, 
that all indexed quantities depend only upon the P populations of neurons (see 
definition ^. ip . i.e. 

9i(t) = f 9p(i)(t) Jij = f Jp{i)p(j) = f I p (i)(t) 

m(t) ~ n p(i) (t) Sj(-) = S pU) (-), (4) 

where x ~ y indicates that the two random variables x and y have the same 
probability distribution. In other words, all neurons in the same population are 
described by identical equations (in law) . 

The g-shapes describe the shape of the postsynaptic potentials and can rea- 
sonably well be approximated by smooth functions. 

In detail we assume that g a , a = 1, • ■ • , P is the Green function of a linear 
differential equation of order fc, i.e. satisfies 

£M*)^r(*) = *(*)■ (5) 

(=0 

where 6(t) is the Dirac delta function. 

The functions bi a (t), I = 0, • ■ ■ , k, a = 1, • ■ ■ , P, arc assumed to be continu- 
ous. We also assume for simplicity that 

6fc Q (i) = c«^0, (6) 

for all t € K, a = 1, ■ ■ ■ ,P. We note D\ the corresponding differential operator: 

k 

dt l 



D k a g Q (t)^Y,b la (t)^(t) = 5(t) (7) 



1=0 



'^When we come to the mean-field equations they will be modeled as random variables. 
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Applying D\ to both sides of ([3]), using ([7]) and the fact that fj(s) = Sj(Vj(s)), 
we obtain a fcth order differential equation for Vi 

N 

DfVi(t) = J2 JijSAVj{t)) + h{t) + m(t). (8) 
j'=i 

With a slight abuse of notation, we split the sum with respect to j into P sums: 

p n p 

0=1 3=1 

We classically turn the fcth-ordcr differential equation ((SJ) into a fc-dimensional 
system of coupled lst-ordcr differential equations (we divided both sides of the 
last equation by c,, see equation ©): 

<*Vk(*) = Vi +li (t)dt l = 0,---,k-2 

(9) 

A well-known example of g-shapes, see section [2.2.31 below or 22j|, is 

ff (t) = tfe-'/^i), (10) 

where K(i) is the Heaviside function. This is an exponentially decaying postsy- 
naptic potential corresponding to 

k=l &i(*) = 4 and 6o(t) = _L 

A AT 

in equation 

Another well-known example is 

ff (t) = Kte-^Yit). (11) 

This is a somewhat smoother function corresponding to 



fc = 2 b 2 (t) = ± h(t) = - b Q {t) = \ 

J\ T T 



in equation (J5]). 



The dynamics 

We modify the equations © by perturbing the first k —1 equations with Brow- 
nian noise and assuming that rij(i) is white noise. This has the effect that the 
quantities that appear in equations © are not anymore the derivatives up to 
order k — 1 of V.. This becomes true again only in the limit where the added 
Brownian noise is null. This may seem artificial at first glance but f) it is a 
technical assumption that is necessary in the proofs of the wcll-poscdness of the 
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mean- field equations, see assumptions 13.11 below, and 2) it generates a rich class 
of external stochastic input, as shown below. With this in mind, the equations 
([5]) now read 



dV H (t) 
dVk-u(t) 



-u(t)dt + f u (t)dWu(t) i = 0, - - - , * - 2 

- Ef=o K(€) (t)Vu(t) + Ej JijS v{j) (v, (t)) + i m (t) 

fk-u(t)dWk-n(t) 



W u (t), I = 0, 



(12) 

k — 1, i = 1, • ■ • , iV, are kN independent standard Brownian 
processes. Because we want the neurons in the same class to be essentially 
identical we also assume that the functions fu(t) that control the amount of 
noise on each derivative satisfy 



/«(*)=/!„«(*), 1 = 0,- 



1. 



1,- 



,N 



Note that in the limit fi a (t) = for I = 0, ■ • • , k — 1 and a = 1, • ■ • , P, 
the components Vu(t) of the vector Vj(£) are the derivatives of the membrane 
potential Vi, for I = 0, • • • , k — 1 and the equations (fT2|) turn into equations ([9|). 
The system of differential equations (fT2|) implies that the class of admissible 
external stochastic input rii{t) to the neuron i are Brownian noise integrated 
through the filter of the synapse, i.e. involving the Ith primitives of the Brownian 
motion for / < k. 

We now introduce the k — 1 A^-dimensional vectors Vj(i) = [Vu, • • • , VJjv] t , 
I = 1, ■ ■ ■ , k — 1 of the Zth-order derivative (in the limit of fi p ^) (t) = 0) of V(i), 
and concatenate them with V(£) into the A^fc-dimensional vector 



V(t) = 



V(t) 
Vi(t) 

V fc _i(t) 



(13) 



The iV-neurons network is described by the iVfc-dimensional vector V(£). By 
definition the Zth A^-dimcnsional component V; of V is equal to V/ . In the limit 
fia{t) = we have 

d l V 
IF 

We next write the equations governing the time variation of the k iV-dimensional 
sub- vectors of V(t), i.e. the derivatives of order 0, ... k — 1 of V(£). These are 
vector versions of equations (fl"2"l) . We write 



V, = V, 



l = 0,---,k-l, with V =V 



dV,(t) = Vj +1 (*) dt + F t (t) ■ dW,(t) I = 0,- ■ ■ , k - 2. 
Fi(t) is the TVxiV diagonal matrix diag(/a(£), ■ • ■ , fa (£),•• • ,fip(t),-- ■ 

S V ' " V- 

JVi iVf 

where fi a (t), a = 1, ■ ■ ■ , P is repeated 7V Q times, and the W;(t), Z = 0, • 
are fc — 1 TV-dimensional independent standard Brownian processes. 



(14) 
,fip(t)), 

■■ ,k-2, 
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The equation governing the (k — l)th differential of the membrane po- 
tential has a linear part determined by the differential operators D^, a = 
1, • ■ ■ , P and accounts for the external inputs (deterministic and stochastic) 
and the activity of the neighbors. We note C(t) the N x Nk matrix describ- 
ing the relation between the neurons membrane potentials and their deriva- 
tives up to the order k — 1 and the (k — l)th derivative of V. This matrix 
is defined as the concatenation of the k N x N diagonal matrixes B;(t) = 
diag(6 a (i),--- ,&u(t),-" ,hp(t),--- ,bip(t)) for I = (),■■■ ,k-l: 



N 1 



C(t) = p (t), 



,B fc _i(t)] 



We have: 

dV fc _i(t) = ( - C(t) ■ V(t) + (J ■ S(V (t))) + I(t)) * 

+ F fc _ 1 (t)-dW fc _i(t) J (15) 

where Wfc_i(f) is an A^-dimcnsional standard Brownian process independent of 
Wi(i), i = 0, ■ ■ ■ , fc— 2. The coordinates of the TV-dimensional vector I(i) are the 
external deterministic input currents, = [Ii(t),--- ,Ii(t),--- ,Ip(t),--- ,Ip(t)] 

^— ' " v ' 

Ni N P 

J the N x N matrix of the weights which are equal to J p (i) P (j) (see equations 
(01), and S is a mapping from ]R N to such that 



S(V)i = S p{i) (Vi) for i = l r ..,N. 



(16) 



We define 



L(t) = 



OATxAf 



B (t) Bi(i) 



®NxN 

Mat 
Bfc_i(t) 



where Id at is the N x N identity matrix and Oatxjv the N x N null matrix. We 
also define the two fciV-dimcnsional vectors: 





Oat 




Oat 




' Oat " 


u t = 


Oat 

. J-5(V (i)) . 




Oat 

. J • S(V(t)) . 


and Ij = 


On 
. I(*) . 



where Oat is the iV-dimensional null vector. 

Combining equations (fl4|) and (JT5J) the full equation satisfied by V can be 
written: 

dV(t) = (-L(i)V(i) + U t + I t J <ft + F(t) ■ dW t , (17) 

where the kN x kN matrix F(t) is equal to diag(Fo, • • • , Fk-i) and W t is an 
fcTV-dimcnsional standard Brownian process. 
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2.2 The Mean-Field equations 



One of the central goals of this paper is to analyze what happens when we let 
the total number N of neurons grow to infinity. Can we "summarize" the kN 
equations (|17[) with a smaller number of equations that would account for the 
populations activity? We show that the answer to this question is yes and that 
the populations activity can indeed be represented by P stochastic differential 
equations of order k. Despite the fact that their solutions are Gaussian pro- 
cesses, these equations turn out to be quite complicated because these processes 
are non-Markovian. 

We assume that the proportions of neurons in each population are non- 
trivial, i.e. : 



If it were not the case the corresponding population would not affect the global 
behavior of the system, would not contribute to the mean-field equation, and 
could be neglected. 

2.2.1 General derivation of the mean-field equation 

When investigating the structure of such mesoscopic neural assemblies as cor- 
tical columns, experimentalists are able to provide the average value Jy of the 
synaptic efficacy Jy of neural population j to population i. These values are 
obviously subject to some uncertainty which can be modeled as Gaussian ran- 
dom variables. We also impose that the distribution of the Jys depends only 
on the population pair a = p(i),/3 = p(j), and on the total number of neurons 
Np of population (3: 



We also make the additional assumption that the Jjy's are independent. 
This is a reasonable assumption as far as modeling cortical columns from ex- 
perimental data is concerned. Indeed, it is already difficult for experimentalists 
to provide the average value of the synaptic strength J a p from population j3 to 
population a and to estimate the corresponding error bars ((J a p), but measuring 
synaptic efficacies correlations in a large assembly of neurons seems currently 
out of reach. Though, it is known that synaptic weights are indeed correlated 
(e.g. via synaptic plasticity mechanisms), these correlations are built by dynam- 
ics via a complex interwoven evolution between neurons and synapses dynamics 
and postulating the form of synaptic weights correlations requires, on theoret- 
ical grounds, a detailed investigation of the whole history of neurons-synapses 
dynamics. 

Let us now discuss the scaling form of the probability distribution (|18|) 
of the Jy 's, namely the division by Np for the mean and variance of the 



lim -f =n a G (0,1) Va G {1, P}, and V 

<!— too /V ^ — ' 



= 1. 



a 




(18) 
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Gaussian distribution. This scaling ensures that the "local interaction field" 
y^fi JijS(Vj(t)), summarizing the effects of the neurons in population /3 on 
neuron i, has a mean and variance which do not depend on Np and is only 
controlled by the phcnomcnological parameters J a f3,<? a f3- 

We are interested in the limit law when N — ► oo of the ./V-dimensional 
vector V defined in equation ^ under the joint law of the connectivities and 
the Brownian motions, which we call the mean- field limit. This law can be 
described by a set of P equations, the mean-field equations. As mentioned in the 
introduction these equations can be derived in several ways, either heuristically 
as in the work of Amari Q, 0] , Sompolinsky [l5|, 34 1 , and Cessac 0, HH , or 



rigorously as in the work of Benarous and Guionnet 0, [H, 0|. We derive 
them here in a pedestrian way, prove that they are well-posed, and provide 
an algorithm for computing their solution. 

The effective description of the network population by population is possible 
because the neurons in each population are interchangeable, i.e. have the same 
probability distribution under the joint law of the multidimensional Brownian 
motion and the connectivity weights. This is the case because of the relations 
(fJJ and (fT6)) which imply the form of equation (fT"7|) . 

The mean ideas of dynamic mean-field equations. 

Before diving into the mathematical developments let us comment briefly what 
are the basic ideas and conclusions of the mean-field approach. Following equa- 
tion ijSJl , the evolution of the membrane potential of some neuron i in population 
a is given by: 

k d l V- N 

E M*)-*r(*) = E JiiSiWV)) + Ht) + rn(t), p(i) = a. (19) 

1=0 3=1 

Using the assumption that Si, Ii, n% depend only on neuron population, this 
gives: 

k i P 

E M*)-^r(*) = E + 7 «(*) + »*(*). 1 e a < ( 2 °) 

1=0 (3=1 

where we have introduced the local interaction field f]i/3(V(t)) — J2^=i JijSpiVjit)), 
summarizing the effects of neurons in population (3 on neuron i and whose prob- 
ability distribution only depends on the pre- and postsynaptic populations a and 
p. 

In the simplest situation where the Jij's have no fluctuations {<J a = 0) this 
field reads rup(V(t)) = J a ^p{V{t)). The term $p(V(t)) - ^ E£i S fj {V 3 {t)) 
is the frequency rate of neurons in population /3, averaged over this population. 
Introducing in the same way the average membrane potential in population f3, 
Vp(t) = ^ Vj(t), one obtains: 
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k i P 

E M*)^r(*) = £ J aj8 */j(V(t)) + / a (t) + n a (t). (21) 
;=o j8=i 

This equation ressembles very much equation (|19|) if one makes the following 
reasoning: "Since &p(V(t) is the frequency rate of neurons in population f3, 
averaged over this population, and since, for one neuron, the frequency rate is 
Vi{t) = Si(Vi(t)) let us write $p(V(t)) = Sp{V {t)f . This leads to: 

fa 1 P 

E M*)^r(*) = E Ja0S (Vp{t)) + I a (t) + n a (t), (22) 

1=0 (3=1 

which has exactly the same form as (fT9|) but at the level of a neuron population. 
Equations such as (|22[) . which are obtained via a very strong assumption: 




are typically those obtained by Jansen-Rit (27j . Surprisingly, they are correct 
and can be rigorously derived, as discussed below, provided a a p = 0. 

However, they cannot remain true, as soon as the synaptic weights fluctuate. 
Indeed, the transition from equation (|T5|) to equation (f2"2"| corresponds to a 
projection from a iVP-dimcnsional space to a P-dimensional one, which holds 
because the NP x NP dimensional synaptic weights matrix has in fact only P 
linearly independent rows. This does not hold anymore if the Jij's are random 
and the synaptic weights matrix has generically full rank. Moreover, the effects 
of the nonlinear dynamics on the synaptic weights variations about their mean, 
is not small even if the o~ a ps are and the real trajectories of (|19p can depart 
strongly from the trajectories of (|2"2"|) . This is the main message of this paper. 

To finish this qualitative description, let us say in a few words what hap- 
pens to the mean-field equations when cr Q| g ^ 0. We show below that the local 
interaction fields T] a p(V(t)) becomes, in the limit Np — * oo, a time dependent 
Gaussian field U a p(t). One of the main results is that this field is non Marko- 
vian, i.e. it integrates the whole history, via the synaptic responses g which 
are convolution products. Despite the fact that the evolution equation for the 
membrane potential averaged over a population writes in a very simple form: 

E M*)-^r(*) - E u <#® + 7 « W + ( 24 ) 

1=0 0=1 

it hides a real difficulty, since U a p(t) depends on the whole past. Therefore, the 
introduction of synaptic weights variability leads to a drastic change in neural 
mass models, as we now develop. 
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The Mean-Field equations We note C{[t , T], R p ) (respectively C((-oo, T], R J 
the set of continuous functions from the real interval [to, T] (respectively (—00, T]) 
to R p . By assigning a probability to subsets of such functions, a continuous 
stochastic process X defines a positive measure of unit mass on C([to, T], R p ) 
(respectively C((— 00, T], R p ))- This set of positive measures of unit mass is 
noted Mf (C*([t ,T],R p )) (respectively M f(C((-oo, T], R p )). 

We now define a process of particular importance for describing the limit 
process: the effective interaction process. 

Definition 2.2 (Effective Interaction Process). Let X G Mf(C([t , T], R p )) 
(resp. VVff (C((— 00, T], R p )) be a given Gaussian stochastic process. The ef- 
fective interaction term is the Gaussian process U x € M^(C([to, T], R PxP )), 
(resp. .M+(C*((-oo,T],R PxP )) defined by: 



E 



U&t) = J a pm x (t) 



<m«^«^ ""V"* " ' <25) 

10 otherwise 
where 

mf(t) d ^ f E[^(^(t))], 

and 

Af(t, S )H f E[^(^(t))^(X^( S )) 



In order to construct the solution of the mean-field equations (see section [3]) 

f(t) and Af 



we will need more explicit expressions for m x (t) and A x (t, s) which we obtain 



in the next proposition. 

Proposition 2.1. Let /j(t) = E [X t ] be the mean of the process X and C(t, s) — 
E \(X t — fi(t))(X s — ^i(s)) T ] be its covariancc matrix. The vectors m x (t) and 
A x (t,s) that appear in the definition of the effective interaction process U x 
are defined by the following expressions: 



m x (t) = J Sp (xjCppfat) + n P {t)j Dx, 



(26) 



and 



A x u x f c f y/Cpp(t,t)Cpp(s,s)-Cpp(t,s) s Cpp^s) 



Sfl (yy/Cpp{t,t)+^(t) \ Dx Dy, (27) 



where 



1 _ai 

Da; = , e 2 dx. 

V2^ 



is the probability density of a zero-mean, unit variance, Gaussian variable. 
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Proof. The results follow immediatly by a change of variable from the fact that 
Xp(t) is a univariate Gaussian random variable of mean np(t) and variance 
C/3p(t, t) and the pair {X p{t) , X p(s)) is bivariate Gaussian random variable with 
mean (fip(t), fxp(s)) and covariance matrix 



Cf3p(t,s) 

Cf}p{s, s) 



□ 



Choose P neurons ii, . . . , ip, one in each population (neuron i a belongs to 
the population a). We define the fcP-dimensional vector V^ N '(t) by choosing, 
in each of the k TV-dimensional components V;(t), I = 0, • • • , k— 1, of the vector 
V(t) defined in equation (|13p the coordinates of indexes £].,••• , ip. Then it 
can be shown, using cither a heuristic argument or large deviations techniques 

(see appendix [A}, that the sequence of fcP-dimensional processes (Vt>to) 

r solution 



converges in law to the process V(t) 
of the following mean-field equation: 



[v(t) T ,W 



,Vn(i)' 



dV(t) 



-L(t) V(t) + + ) dt + F(t) ■ dW 



(28) 



L is the Pk x Pfc matrix 



L(<) = 



OpxP Idp 

OpxP OpxP 

B (t) Bx(t) 



OpxP 
OpxP 

Idp 

B fc _i(t) 



The P x P matrixes Bj(f), Z = 0, • • • , — 1 are, with a slight abuse of notations, 
equal to diag(&/i(£), • • • , bip(t)). (W t ) t > to is a fcP-dimensional standard Brow- 
nian process. U v has the law of the P-dimensional effective interaction vector 
associated to the vector V (first P-dimensional component of V) and is statis- 
tically independent of the external noise (W t ) 4 > to and of the initial condition 
V(fo) (when to > — oo): 



Op 




' Op " 




I(t) = 




Op 




Op 


_ ■ 1 _ 




. m _ 



F(t) = diag(F (i), • • • .PWi)). 



Wc have used for the matrixes Fi(t), I — 0. ■ • ■ ,k- 
as for the matrixes Bj(i), i.e. f = diag(/n(f), ■ 
I(t) is the P-dimensional external current [ii(£) • 



1 the same abuse of notations 
•• ,/ IP (i))forJ = 0,-.. ,fe-l. 
■Jp(t)] T . 



The process (U^) t > to is a P x P-dimensional process and is applied, as a 



matrix, to the P-dimensional vector 1 with all coordinates equal to 1, resulting 
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in the P-dimensional vector • 1 whose mean and covariance function can be 
readily obtained from definition 12.21 

p 

E[uM]=X)J a/S n#(i), m}(t)=E[S (V (t))] (29) 

13=1 

and 

f EjU^/J^CM) if a = 7 
Cov((Ur • l) a (U]f ■ 1) 7 ) = ^ (30) 

[ otherwise 

We have of course 

AV(t >a )=E[5 /J (V /3 (t))5 /3 (V^( a ))] 

Equations (|28[) are formally very similar to equations (|17[) but there are some 
very important differences. The first ones are of dimension kP whereas the 
second are of dimension kN which grows arbitrarily large when N — ► oo. The 
interaction term of the second, J • S(V(t)), is simply the synaptic weight matrix 
applied to the activities of the N neurons at time t. The interaction term of 
the first equation, U^, though inocuous looking, is in fact quite complex (see 
equations (|2"9")) and (|3"0"1) 1. In fact the stochastic process Uj , putative solution 
of equations (|28|) , is in general non Markovian. 

To proceed further we formally integrate the equation using the flow, or 
resolvent, of the equation (|2"5|) . noted $L(t, to) (see appendix[B]), and we obtain, 
since we assumed L continuous, an implicit representation of V(t): 



V(t) = M*.*o)V(*b)+ / • (Uy+I(a) )ds 



f $ L (t,s)-F 



(s) ■ dW s (31) 



We now introduce for future reference a simpler model which is quite fre- 
quently used in the description on neural networks and has been formally an- 
alyzed by Sompolinski and colleagues, [TBI, 34 1 in the case of one population 
(P=l). 



2.2.2 Example I: The Simple Model 

In the Simple Model, each neuron membrane potential decreases exponentially 
to its rest value if it receives no input, with a time constant r a depending only on 
the population. In other words, we assume that the g-shape describing the shape 
of the PSPs is equation (TTU)) , with K = 1 for simplicity. The noise is modeled 
by an independent Brownian process per neuron whose standard deviation is 
the same for all neurons belonging to a given population. 
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Hence the dynamics of a given neuron i from population a of the network 
reads: 



dVi(t) 



Vj(t) 



P 



+ j ^ s pu) (W*)) + ^(o (*) 



/3=lj=l 



/ pW dWi(t). (32) 



This is a special case of equation (|12|l where fc = 1, bo a (t) = l/r , &i a (t) = 1 
for a = 1, • • • , P. The corresponding mean-field equation reads: 



p 

dV a (t) = f-Y^l+J2 U^(t)+I a (t)) dt+f a dW a (t), Va e {1, . . . , P}, (33) 

Tq /3=1 

where the processes (W a (t)) t >t are independent standard Brownian motions, 
U v (£) = (U^g(t); a,fl £ {1, . . . , P}) t is the effective interaction term, see 
definition ^. 21 This is a special case of equation l|28p with L = diag(^-, • • • , i), 
and F = diag(/i,- • • , /». 

Taking the expected value of both sides of equation (|33|) and using we obtain 
equation (|26f that the mean ^ Q (t) of V a (t) satisfies the differential equation 



dn a (t) fi a (t) 



dt Tr, 

13=1 



J2 J ^ J ^0 (xyfc^(t~£) + ftfi(t)j Dx + I a (t), 



If Cpp(t,i) vanishes for all t > to this equation reduces to: 



djj, a {t) fi a (t) 



p 



^2 J *P S fi + Ia{t), 



dt T a 

/3=1 

which is precisely the "naive" mean-field equation (|22|) obtained with the as- 
sumption (j23|) . We see that equations (|22|) are indeed correct, provided that 

Cft8(t,t) = 0, Vt > t . 

Equation ()33[) can be formally integrated implicitly and we obtain the fol- 
lowing integral representation of the process V a (t) : 



V a (<) = e^'"'^* V«(*o) + f e- ( *- s)/r ° ( ^ C7^(s) + I a {s)) ds 

+ fa f e- (t - s)/T "dW a (s) (34) 

Jto 

where to is the initial time. It is an implicit equation on the probability distribu- 
tion of V(i), a special case of (01]), with $ L (t, t ) = diag(e- ( *- to)/ri , • • • , e -(*-* o) / Tp ). 
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The variance C aa (t,t) of V a (t) can easily be obtained from equation (f3~T| . 
It reads 



C aa (t,t) = er 2t /^ [c aa (t ,t ) + ^ (e- - l) 



3=1 



ctL / / e (u+v)/r «A p (u,v)dudv 

to J to 

where A ( a(it,w) is given by equation (1271) . 

If cr Q /3 = and if s a = then C aa (t,t) = 0,Vi > io is a solution of this 
equation. Thus, mean-field equations for the simple model reduce to the naive 
mean-field equations (|22[) in this case. This conclusion extends as well to all 
models of synaptic responses, ruled by equation j5]). 

However, the equation of C aa (t, t) shows that, in the general case, in order 
to solve the differential equation for /i Q (£), we need to know the whole past of 
the process V. This examplifics a previous statement on the non Markovian 
nature of the solution of the mean-field equations. 



2.2.3 Example II: The model of Jansen and Rit 

One of the motivations of this study is to characterize the global behavior of an 
assembly of neurons in particular to get a better understanding of recordings of 
cortical signals like EEG or MEG. One of the classical models of neural masses 
is Jansen and Rit's mass model [27|, in short the JR model (see figure [2|). 

The model features a population of pyramidal neurons that receives in- 
hibitory inputs from local inter-neurons, excitatory feedbacks, and excitatory 
inputs from neighboring cortical units and sub-cortical structures such as the 
thalamus. The excitatory input is represented by an external firing rate that has 
a deterministic part I\ (t) accounting for specific activity of other cortical units 
and a stochastic part ni(t) accounting for a non specific background activity. 
We formally consider that the excitatory feedback of the pyramidal neurons is 
a new neural population, making the number P of populations equal to 3. We 
also represent the external inputs to the other two populations by the sum of a 
deterministic part Ij(t) and a stochastic part nj(t), j = 2, 3, see figure [2] 

In the model introduced originally by Jansen and Rit, the connectivity 
weights were assumed to be constant, i.e. equal to their mean value. Nev- 
ertheless, there exists a variability of these coefficients, and as we show in the 
sequel, the effect of the connectivity variability impacts the solution at the level 
of the neural mass. Statistical properties of the connectivities have been studied 
in details for instance in Q. 

We consider a network of ./V neurons, N a , a = 1, 2, 3 belonging to population 
a. We index by 1 (respectively 2, and 3) the pyramidal (respectively excitatory 
feedback, inhibitory interneuron) populations. We choose in each population 
a particular neuron indexed by i a , a = 1,2,3. The evolution equations of the 
network can be written for instance in terms of the potentials , Vi 2 and Vi 3 
labelled in figure [2] and these equations read: 
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(a) Populations involved in Jansen's model (b) Block diagram 



Figure 2: a. Neural mass model: a population of pyramidal cells interacts 
with itself in an excitatory mode and with an inhibitory population of inter- 
neurons, b. Block representation of the model. The g boxes account for the 
synaptic integration between neuronal populations. S boxes simulate cell bodies 
of neurons by transforming the membrane potential of a population into an 
output firing rate. The coefficients J a p are the random synaptic efficiency of 
population f3 on population a (1 represents the pyramidal population. 2 the 
excitatory feedback, and 3 the inhibitory inter- neurons) . 



N 2 N 3 

= 5i * I £ JiiiS(Vj) + J2 J idS(Vj) + h 4 

JVi 

= 9i * | J ^jS(Yi) + h + n 2 

JVi 

Vis =93* J ^J S ( V l) + h + "3 

In the mean- field limit, denoting by V Q , a = 1, 2, 3 the average membrance 
potential of each class, we obtain the following equations: 



Vi = gi * (U 
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u 



V 2 = gi * (f7 21 + J 2 + n 2 ) 
V 3 = g 3 * (C/ 3 V i + h + n 3 ) 



(35) 



where U v = (C/2/3)a,/3=i-2,3 is the effective interaction process associated with 
this problem, i.e. a Gaussian process of mean: 
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'E[[7}y =Ji a E[5(Va)] 

Ek V 3 ] =Ji 3 E[5(V 3 )] 

E[E7&] =J 21 E[5(Vi)] 

[E[C/V] =J 31 E[5(Vi)], 

All other means correspond to the non-interacting populations and are equal to 
zero. The covariance matrix can be deduced from 

Cov(t/^(i), U^s)) 

where 



ff*aA£(*,«) if a = and 1 = 5 
otherwise 



A^(t, a ) = E[5(V a (*))S(V / ,(«))] 

This model is a voltage-based model in the sense of Ermentrout [l9[ ■ Let us now 
instantiate the synaptic dynamics and compare the mean-field equations with 
Jansen's population equation^ (sometimes improperly called also mean-field 
equations) . 

The simplest model of synaptic integration is a first-order integration, which 
yields exponentially decaying post-synaptic potentials: 



9(t) 



Ke~~ t > 
t < 



Note that this is exactly equation (flQ|) . The corresponding g-shape satisfies the 
following lst-order differential equation 

g(t)=--g(t) + KS(t), 

T 

In this equation r is the time constant of the synaptic integration and K the 
synaptic efficiency. The coefficients K and r are the same for the pyramidal and 
the excitatory feedback population (characteristic of the pyramidal neurons and 
defining the g-shape g\), and different for the inhibitory population (defining 
the g-shape gs). In the pyramidal or excitatory (respectively the inhibitory) 
case we have K = K±, r = T\ (respectively K = K 3 , t = t 3 ). Finally, the 
sigmoid functions S is given by 

S(v) 



1 + e r(v -v) 



where v max is the maximum firing rate, and vq is a voltage reference. 

With this synaptic dynamics we obtain the first-order Jansen and Rit's equa- 
tion: 

• = _ J_ Vl + K x (U? 2 + U? 3 + ni ) 
d£ = -±V 3 + K 1 (UZ l +I a +n 2 ) . (36) 



^ =-±V 3 + K 3 (uV+l 



3 + n 3 ) 



J We have modified the original model which is not voltage-based. 



20 



The "original" Jansen and Rit's equation [23j, |27j amount considering only the 
mean of the process V and assuming that E[Si(Vj)] = 5j(E[Vj]) for i,j £ 
{1,2,3}, i.e. that the expectation commutes with the sigmoidal function S. 
This is a very strong assumption, and that the fluctuations of the solutions 
of the mean-field equation around the mean imply that the sigmoid cannot be 
considered as linear in the general case. 



A higher order model was introduced by van Rotterdam and colleagues 37 
to better account for the synaptic integration and to better reproduce the char- 
acteristics of real postsynaptic potentials. In this model the g-shapes satisfy a 
second order differential equation: 



Kte~- t > 
t < 



<rvi 


_ _2_ 


dVi 


dt 2 


Tl 


dt 


d 2 V 2 


_ _2_ 


dV 2 


dt 2 


Tl 


dt 


d 2 V 3 


_ __2_ 


dV 3 


dt 2 


T3 


dt 



We recognize the g-shapc defined by equation (|11[) solution of the second-order 
differential equation y(t) + -y(t) + -^y{t) = K8{t). With this type of synaptic 
integration, we obtain the following mean-field equations: 

- ^Vi + K x (U? 2 + U& + h+ m) 
--i v 2 + ^i(C/2 V i+^2 + n 2 ) (37) 
-^V 3 + K 3 (uV+I 3 + n 3 ) 

Here again, going from the mean-field equations (|37|) to the original Jansen 
and Rit's neural mass model consists in studying the equation of the mean 
of the process given by (|37|) and commuting the sigmoidal function with the 
expectation. 

Note that the introduction of higher order synaptic integrations results in 
richer behaviors. For instance, Grimbert and Faugeras [23} showed that some 
bifurcations can appear in the second-order JR model giving rise to epileptic 
like oscillations and alpha activity, that do not appear in the first order model. 



3 Existence and uniqueness of solutions in finite 
time 

The mean-field equation (|31| is an implicit equation of the stochastic process 
(V(t))t>t - We prove in this section that under some mild assumptions this 
implicit equation has a unique solution. These assumptions are the following. 

Assumption 3.1. 



1. The matrix L(t) is C° and satisfies ||L(t)|| < for all t in [io,? 1 ], for 
some matrix norm || j| and some strictly positive constant k^. 
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2. The matrix F(t) has all its singular values lowerbounded (respectively 
upperbounded) by the strictly positive constant A^j n (respectively Aj^ax) 
for all t in [t ,T]. 

3. The deterministic external input vector l(t) is bounded and we have 
PWIloo < /max for ah t in [io^P] and some strictly positive constant 

/max 

This solution is the fixed point in the set M.f(C([to, T], lR fep )) of fcP-dimcnsional 
processes of an equation that we will define from the mean-field equations. We 
will construct a sequence of Gaussian processes and prove that it converges in 
distribution toward this fixed point. 

We first recall some results on the convergence of random variables and 
stochastic processes. 

3.1 Convergence of Gaussian processes 

We recall the following result from 0] which formalizes the intuition that a 
sequence of Gaussian processes converges toward a Gaussian process if and 
only if the means and covariance functions converge. In fact in order for this 
to be true, it is only necessary to add one more condition, namely that the 
corresponding sequence of measures (elements of Mf(C([to,T],]R, kP ))) do not 
have "any mass at infinity" . This property is called uniform tightness [6(] . More 
precisely we have 

Definition 3.1 (Uniform tightness). Let {X n }^ =1 be a sequence of fcP-dimensional 
processes defined on [to, T] and P n be the associated elements of M.f (C([to, T], R fcp )). 
The sequence M.^ (C([fo, T], R fep )) is called uniformly tight if and only if for all 
e > there exists a compact set K of C([to,T], R fep ) such that P n (K) > 1 — e, 
n > 1. 

Theorem 3.1. Let {X n \ < ^L l be a sequence of kP- dimensional Gaussian pro- 
cesses defined on [to,T] or on an unbounded interval of Hi. The sequence con- 
verges to a Gaussian process X if and only if the following three conditions are 
satisfied: 

• The sequence {X n \ < ^L l is uniformly tight. 

• The sequence /J. n (t) of the mean functions converges for the uniform norm. 

• The sequence C' n of the covariance operators converges for the uniform 
norm. 

We now, as advertised, define such a sequence of Gaussian processes. 

Let us fix Zo, a fcP-dimensional Gaussian random variable, independent of 
the Brownian and of the process ((X) t )t^[t 0y T]- 

3 We note T(t) the matrix F(t)F(t) T . 

4 In 0, Chapter 3.8], the property is stated whenever the mean and covariance are defined 
on a separable Hilbert space. 
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Definition 3.2. Let X be an element of M\ (C([t , T], R fcp )) and T k be the 
function Mf(C([t ,T},M kP )) -> 7M+(C([i , T],~R kP )) such that 

JF fe (X) t = $ L (Mo) • + / • (Uf +I(s)) 

Jt a 

+ f $ L (t,s)-F(s)dW s 

J t 

where and I(s) are defined^ in section [5] 

Note that, by definition 12.21 the random process (Tk(X)) tG [ tg x], k > 1 is 
the sum of a deterministic function (defined by the external current) and three 
independent random processes defined by Zq, the interaction between neurons, 
and the external noise. These three processes being Gaussian processes, so is 
(^k{^))te[t ,T]- Also note that Tk{X) t<) = Zq. It should be clear that a solution 
V of the mean-field equation (f3"Tj) satisfies V(to) = Zq and is a fixed point of Tk, 
i.e. T k (V) t = V(t). 

Let X be a given stochastic process of Mi(C([to, T], R fcp )) such that X to = Z 
(hence X to is independent of the Brownian). We define the sequence of Gaussian 
processes {X n }™ =0 e Mf(C([t , T], TR kp ) by: 

= X 

= n(x n ) = ri n \x ). n>o, 4 0) = u (38) 

In the remaining of this section we show that the sequence of processes 
jjrWpQ joo_ Q converges in distribution toward the unique fixed-point Y of Tk 
which is also the unique solution of the mean-field equation (|3"Tj) . 




3.2 Existence and uniqueness of a solution for the mean- 
field equations 

The following upper and lower bounds are used in the sequel. 

Lemma 3.2. Consider the Gaussian process ((Uf- • l)t)te[t ,T] ■ U A is defined 
in 12.21 and 1 is the P-dimensional vector with all coordinates equal to 1 . We 
have 

||E [Uf • 1] L < \i = max]T | J Q/3 | U^IU (39) 



for all to < t < T. The maximum eigenvalue of its covariance matrix is up- 
perbounded by o^ lax = f max Q J2a a \a Halloo where ||S/3||oo is the supremum of 
the absolute value of Sp. We also note er^ in = f min aj/ 3 cr^. 

Proof. The proof is straightforward from definition 13.21 □ 

5 For simplicity we abuse notations and identify X and X. 
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The proof of existence and uniqueness of solution, and of the convergence 
of the sequence (|38)) is in two main steps. We first prove that the sequence of 
Gaussian processes {Tl (X)}%L , > 1 is uniformly tight by proving that it 
satisfies Kolmogorov's criterion for tightness. This takes care of condition 1) 
in theorem 13.11 We then prove that the sequences of the mean functions and 
covariancc operators are Cauchy sequences for the uniform norms, taking care 
of conditions 2) and 3). 



3.2.1 Uniform tightness 

We first recall the following theorem due to Kolmogorov [i^, Chapter 4.1]. 

Theorem 3.3 (Kolmogorov's criterion for thightness). Let {X n }^ =1 be a se- 
quence of kP -dimensional processes defined on [to,Tj. If there exist a, (3, C > 
such that 



E 



\\X n (t) - X n (s)f <C\t~s\ 1+a Vs,t€[t ,T) n>l, 



then the sequence is uniformly tight. 

Using this theorem we prove that the sequence {^^(X)}^^, k > 1 sat- 
isfies Kolmogorov's criterion for (3 = 4 and a > 1. The reason for choosing 
= 4 is that, heuristically, dW ~ (dt) 1 / 2 . Therefore in order to upperbound 

E — X n (s)||' 3 by a power of \t — ,s| greater than or equal to 2 (hence 

strictly larger than 1) we need to raise ||X„(t) — X n (s)|| to a power at least 
equal to 4. The proof itself is technical and uses standard inequalities (Cauchy- 
Schwarz' and Jensen's), properties of Gaussian integrals, elementary properties 
of the stochastic integral, and lemma 13.21 It also uses the fact that the input 
current is bounded, i.e. that max a= i.... .psup t g[j o T ] | Jq, (t) | < I max < oo, this is 
assumption 3 in 13.11 



Theorem 3.4. The sequence of processes < (X) > , k > 1 is uniformly 

I J n— 

tight. 

Proof. We do the proof for k = 1, the case k > 1 is similar. If we assume that 
n > 1 and s < t we can rewrite the difference J r [ n \x) t — J r [ n \x) s as follows, 
using property i) in proposition IB. II in appendix [Bj 

r[ n \x) t - T[ n \x) s = ($i(t,t ) - $ L (s,t ))X to 

+ (<S> L (t,s)-Id) ^Us^^u 1 ~ {X) -ldu+ Ql^^Uu 1 {X) -ldu 

J t J s 

+ ($ L (t,s)-Id) f $ L (s,u)F(u)dW u + f <S> L (t,u)F(u)dW u 

Jto Js 

+ ($ L (t,s) - Id) I $ L (s,u)I(u)du+ I $ L (t,u)I(u)du 

J in J S 
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The righthand side is the sum of seven terms and therefore (Cauchy-Schwarz 
inequality): 



\\r[ n \x) t - ^rw.ir < \\*L(t,t ) - <s> L (s,t )r \\x tB \\ 



r 



+ (s-t )\\<£ L (t, S )-Id\\ 2 / \\^ L (s,u)\\ 2 \\\J^"~ {X) -lfdu 

J t 

ft 



(t-s) I \\^ L (t,u)\\ 2 \\U^ 1>iX) -lfdu 



\\$ L (t,s)-ld\\ 



<S> L (s,u)F{u)dW, 



2 




+ 





rt 2 

u)F(u) dW u 



(s-to) 2 ||*L(t,s)-M|| 2 4 a x sup \\$l(s,u) 

ue [t , s] 



ue[s,t] 

Because ||^z-(*, *o) — $l( s ,*o)II < \t — s|||Lj| we see that all terms in the 
righthand side of the inequality but the second one involving the Brownian 
motion are of the order of (t — s) 2 . We raise again both sides to the second 
power, use the Cauchy-Schwarz inequality, and take the expected value: 



1 



E 



\\J?\x)t-Ji n) (X)st < ||$l(Mo)-$l(Mo)|| 4 E[||XJ| 4 ] 



+ (s- t ) 3 \\<S> L (t 7 s) - Jd|| 4 / \\$ L (s, u)\\ 4 E 
ft 

\3 I ii^ n ..mi4 



iu? n " 1)w -i||* 



du 



+ {t-s) A / ||$ L (i, U )|| 4 E 



U?" " W • ill 4 



\4> L (t,s)-Id\\ 4 E 



$ L ( S ,«)F(«)dW l 



du 
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$ L (t,u)F{u)dW u 



+ (s~t ) 4 \\^ L (t,s)-Id\\ 4 lL, sup \\<t> L (s, u)\\ 4 

«e[to,s] 



(40) 



+ (t-^lL x sup \\$ L (t,u) 



Remember that U„ 1 • 1 is a P-dimensional diagonal Gaussian process, 

noted Y u in the sequel, therefore: 

E [ii Y "ii 4 ] = E E [ y » 4 ] + E E K»] E [^ 2 »] • 
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The second order moments are upperbounded by some regular function of fi 
and cr max (denned in lemma 13 - 2[) and, because of the properties of Gaussian 
integrals, so are the fourth order moments. 



We now define B(u) = $i(s, u)F(u) and evaluate E 
We have 



/•B(u)dW„ 



E 



B{u)dW u 



= E 



B(u)dW u 



In 



E 



p p 



E E B l3 { U )dwi e / B ^) dW u 
i=i \j=i Jt ° j \k=i Jt ° 

e e [ f («) ^ r B tlkl («) ^ r s i2J2 («) dw» r s l2fc2 ^ 

Because is by construction independent of W£ Ifi^j and E f* Bij (u)dWl 

for all i, j (property of the It 6 integral), the last term is the sum of only three 
types of terms: 

1 . If ji = ki = j 2 = fc 2 we define 

B HJ (u)dWi) (( B i2j (u)dWl 

to J \Jt 



i : y2 • 

and, using Cauchy-Schwarz 

Ti < E E 

ii,*2 

2. If ji = fci and j'2 = &2 but i ^ j 2 we define 
T 2 = Yl E ([ B ilh {u)dWt 





1/2 






E 











1/2 



which is equal, because of the independence of W^ 1 and W^ 2 to 



E 

il,l2,jl=£j2 



B ilh (u)dWV 



B i2h {u)dWt 



t 



3. Finally, if ji = j 2 and fci = k 2 but ji 7^ k\ we define 



T 3 = E E 



B^uJdW* / B^uJdW* / B llkl (u)dW^ / S i2fel ( u )dW„ fcl 

to J to J to 
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which is equal, because of the independence of Wl 1 and W^ 1 to 



3 



B hh {u)dWt / B. l2n (u)dW^ 



E 



B nkl (u)dW^ / B i2kl ( U )dW^ 



Because of the properties of the stochastic integral, J t " B il j (u) dWl = 7V(0, y J* t * Bf^ (u) duj ) 
hence, because of the properties of the Gaussian integrals 



V2, 



E 



B hj {u)dWi 



to 



for some positive constant k. This takes care of the terms of the form T\. Next 
we have 



E 



B tlJ1 (u)dW^ 



B \n (u) du, 



which takes care of the terms of the form T2- Finally we have, because of the 
properties of the ltd integral 



E 



to 



B^{u)dWt / B l23l (u)dW^ 



to 



B lljl (u)Bi 2jl (u)du, 



which takes care of the terms of the form T*. 



This shows that the term . 



f s <5> L (t,u)F{u) dW u 



in (l40l) is of the order 



of (t — s) 1+a where a > 1. Therefore we have 

E \\\T[ n) (X) t - T < L n \x) s \\ 4 ] < C\t - s| 1+a , a > 1 

for all s, t in [£o,T], where C is a constant independent of t, s. According to 
Kolmogorov criterion for tightness, the sequence of processes \ (X) > is 

I J n— 

uniformly tight. 

The proof for .7-^, fc > 1 is similar. □ 

3.2.2 The mean and covariance sequences are Cauchy sequences 

Let us note /J, n (t) (respectively C n (t,s)) the mean (respectively the covariance 
matrix) function of X n = !F k {X n _i), n > 1. We have: 



/i"(£) = <J>l(Mo)m Zo 
/ 



f <t> L (t,u) (E 




+ !(«)) 


'to V 















n T ... n T 

Up, ) Up, 






+ I(tt) du, (41) 




_ 


«=!,■•• ,P- 
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where m^ n (it) is given by equation (|26p . Similarly we have 

/•Ms 

C n+1 (t,s) = $ L (t,i )£ Zo $ L (Mo) T + / $l(^)F( U )F( U ) t $ l (s, U ) t c^+ 
/ /"^i^uJCovfu^.U^^^wfdudt; (42) 

J t Jto V 7 

Note that the kP x kP covariance matrix Cov yU^" , U^ n J has only one nonzero 
P x P block: 

Cov (u*« , ) ^ = Cov (U*« ■ 1, Uf » ■ 1) , (43) 
According to definition 12.21 we have 




where A*™ (u, v) is given by equation (|2"7|) and Dx is defined in proposition ^. 11 
In order to prove our main result, that the two sequences of functions (fi n ) 
and (C n ) are uniformly convergent, we require the following four lemmas that 
we state without proofs, the proofs being found in appendixes [El to [HI The first 
lemma gives a uniform (i.e. independent of n > 2 and a = 1, • • • , kP) strictly 
positive lowerbound for C™ a (t,i). In what follows we use the following nota- 
tion: Let C be a symmetric positive definite matrix, we note A^ in its smallest 
eigenvalue. 

Lemma 3.5. The following upppcrbounds are valid for all n > 1 and all s, t £ 
[to,T]. 

||M n WIL < ^ (T " t0) m [2o]IL + (M + W)(T - ^)] = /W, 

HC"(*^)lloo <e^ +k ^ T ~^ [p(Z z °) + \l iax (T-t ) + a^ x (T-to) 2 ] ^ S max , 
where and cr max are defined in lemma [3~2[ Aj^ ax is defined in 13. II 

Lemma 3.6. For all t £ [to, T] all a = 1, ■ • ■ , £P, and n > 1, we have 

where A m i n is the smallest singular value of the positive symmetric definite 
matrix <I> L (i, t )& L (t 7 t ) T for t £ [to,T] and A^ s ° is the smallest eigenvalue of 
the positive symmetric definite covariance matrix T, z ° . 

The second lemma also gives a uniform lowerbound for the expression C™ a (s, s)C' t 
C™ a (t,s) 2 which appears in the definition of C n+1 through equations f43|) and 
((27)) . The crucial point is that this function is 0(\t — s|) which is central in the 
proof of the third lemma 13.81 



28 



Lemma 3.7. For all a = 1, • • • , kP and n > 1 the quantity C™ a (s, s)C™ a (t, t) — 
C™ a (t,s) 2 is lowerbounded by the positive symmetric function: 

ni ,\ def I , _ i,2 \S Z ° \T 

where is the strictly positive lower bound, introduced in l3.ll on the singular 
values of the matrix F(u) for u E [to, T]. 

The third lemma shows that an integral that appears in the proof of the uni- 
form convergence of the sequences of functions (p n ) and (C ra ) is uppcrboundcd 
by the nth term of a convergent series. 

Lemma 3.8. The 2n-dimcnsional integral 



pi{ui,vi)\ \ " / Pn-l(Un-l,V n -l) 

[t ,tWs] 2 \ J[to,«lVt.!] 2 \J[to,«„-2V«„_ 2 ] 2 

Pn(u n ,v n )du n dv n \ du n _idv n -i ••• \duidvi, 

'[t ,«„_ 1 Vl)„_ 1 ] 2 J J J 

where the functions pi (v,i , Vi) , i = 1 , • • • , n are either equal to 1 or to 1 / y/ 9{ui, m) 
(the function 6 is defined in lemma [3~7]) . is uppcrboundcd by k n /(n — 1)! for 
some positive constant k. 

With these lemmas in hand we prove proposition 13. 91 The proof is technical 
but its idea is very simple. We find upperbounds for the matrix infinite norm 
of C*" +1 (M) - C n (t,s) and the infinite norm of p n+1 {t) - p n {t) by applying 
the mean value theorem and lemmas 13.61 and 13.71 to the these norms. These 
upperbounds involve integrals of the infinite norms of C n (t, s) — C ,l_1 (i, s) and 
p, n (t)~ and, through lemma [3~71 one over the square root of the function 

9. Proceeding recursively and using lcmma [3T5I one easily shows that the infinite 
norms of C" +1 — C n and p n+1 — p n arc uppcrboundcd by the nth term of a 
convergent series from which it follows that the two sequences of functions are 
Cauchy sequences, hence convergent. 

Proposition 3.9. The sequences of covariance matrix functions C ra (i, s) and 
of mean functions p n (t), s, t in [ioj^l are Cauchy sequences for the uniform 
norms. 

Proof. We have 

C n+1 (t , s)-C n (t,s) = j* £* L {t,u) (Cov (U*« , Uf ») - 

Cov^U^- 1 ,!!^"- 1 ) ^ L {s,v) T dudv. 

We take the infinite matrix norm of both sides of this equality and use the upper- 
bounds H^lMIL < ell L IU( T -*°) = k L and ||$ L (t,w) T || oo < e ll LT IL.( T - to ) = 
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k L T (see appendix |B|) to obtain^ 



\C n+1 (t,s)-C n (t,s)\\ <k L k 



to J to 



Cov(U£»,U?») -Cov(U^-,U^ 



du dv 



kr,k 



LK L T 



to J to 



Cov (U?« • 1, U*« ■ 1) - Cov (U^- 1 • 1, - 1 • 1) 



c?w dv. 



(44) 



According to equations ([27]) we are led to consider the difference j4 n — A n -i, 
where: 



A, 



del' 



S/J ( *+ , V+f$(v)) S (yJc^u) + ^(u) 

d ^ f S [P$(u,v)x+S%(u,v)y + fi(v)] S p [T p n (u)y + fi(uj\ . 

We write next: 

K - A n -i = S [P$(u, v)x + S$(u, v)y + ^(v)] 

(Sp [7?(u)y + i$(u)] S P \T£-\u)y + ^(u)] ) + 
Sp[T^-\u)y + ^f\u) 
(s [P%(u,v)x + S%(u,v)y + i$(v)}-S [v^ 1 (u, v)x + S%~ 1 (u, v)y + ^ («)]) . 
The mean value theorem yields: 

I A„ - A n _x |< H^IL HS^IL ( | x 1 1 V$(u,v) V n f\u, v) | + 
\y\\ S%(u,v)-S%- 1 (u,v) | + | ^(v)-» n f\v) \ + \y\\ VM-^'M I + 

\^(u)-n}-\u)\ 



6 The notation || ||" is introduced in appendix [C] 
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Using the fact that J_ oa \ x \ Dx = y we obtain: 



\C n+1 (t,s)-C n (t,s)\\ oo <k L k LT k c \ 

ft rs 



V n {u,v) ~V n -\u,v)\\ dudv+ 



to J to 



S n {u,v)-S n - L (u,v)\\ dudv+ 



to J to 



(t-to) / ||M»-M n ">)L <*«+(* -to) / -/x"- 1 ( W )|| oo ^+ 

Jt J t 



-{a -to) T n (u)-T n -\u) \ du , 



where the constants fci and /c^t are defined in appendix [Bl and 



(45) 



A similar process applied to the mean values yields: 

^"(uJ-^-^uJll^du), 



where is defined in lemma [3~^1 We now use the mean value theorem and lem- 
mas l3.6l and l3.7l to find upperbounds for ||p n (u, i;) — V n ~ 1 (u 7 w) , ||5™(u,w) — 5™ _1 (u,w)| 
and 1 1 T" (u) - T" _1 (w)|| . We have 



\T?{u)-T;r\u 



where fco is defined in lemma 1X61 Hence: 

1 



< 



1 



2Jkr] 



C n g0 {u,u)-C^\u,u) 



T» - T n ~ L (u)\\ < —= C n (u,u) - C"- A (u,u) 



2Vko 

Along the same lines we can show easily that: 

\\S n (u,v) - 5"- 1 K U )|| oo < fc( \\C n (u,v) - C^fotOlL + 



and that: 



\V n (u,v)-V n ~\u,v)\\ oo <^=^=(\\C n (u,v)-C n - 1 (u,v)\\ oo + 

\\C n (u, u) - C"" 1 ^, u)^ + \\C n (v, v) ~ C n ~\v, v)^ ) , 
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where 9(u,v) is defined in lemma [3~71 Grouping terms together and using the 
fact that all integrated functions are positive, we write: 



\C n+1 (t,s)-C n (t, s )\\ oc < 
1 



to,*Vs] 2 \/0(u,v) 



C n (u,v)-C n - L (u,v)\\ dudv+ 



"C n (u,u) - C*"" 1 ^,?/)!! dudv- 



[t ,tvs] 2 \/0(u 7 v) 

\\C n (u,v) - C"-- 1 ^^)^ dudv+ 

[t ,*Vs] 2 

||C n (u,u)-C ,n - 1 (u,iO|| oo dudv+ 

[to,*Vs] 2 

/ \\f{u)-ii n -Hu)\\ x dudv). (46) 

Note that, because of lemma l3~6l all integrals are well-defined. Regarding the 
mean functions, we write: 



l^n+i^.^H < A / \\C n (u > u)-C n - 1 (u,u)\\ oo dudv+ 

\ J[t ,tVs] 2 

W^W-^-^u^dudv). (47) 

[*o,tVs] 2 / 

Proceeding recursively until we reach C° and fj, we obtain an upperbound 
for \\C n+1 {t,s) - C n (t,s) (respectively for \\n n+1 (t) - /**"(*) H^) which is the 
sum of less than 5" terms each one being the product of k raised to a power 
less than or equal to n, times 2/^ max or 2E max (upperbounds for the norms 
of the mean vector and the covariance matrix defined in lemma [3~5]) . times a 
2n-dimcnsional integral /„ given by 



/ ' ' / p„_i(u„_i,w n _i) 

[t ,tVs] 2 \ Jlto.UiV-Ui] 2 \-'[t0,«r 1 -2V'U„_ 2 ] 2 

Pn{u n ,v n )du n dv n du„_idw„_i ••■ \du\dv\, 

[t ,«„-lVl,„_ 1 ] 2 / / / 



where the functions pi (it,, Vi) , i = 1, • ■ ■ , n are either equal to 1 or to 1 / \J 6{y,i, vi) 
According to lemma 13. 8[ this integral is of the order of some positive con- 
stant raised to the power n divided by (n — 1)!. Hence the sum is less than 
some positive constant k raised to the power n divided by (n — 1)!. By taking 
the supremum with respect to t and s in [ioi^" 1 ] we obtain the same result for 
||C" +1 — C n || (respectively for — jit n || ). Since the series X)n>i ^ s 
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convergent, this implies that \\C n+p — C"!^ (respectively ||^ n+p — ^" I loo) can 
be made arbitrarily small for large n and p and the sequence C n (respectively 
H n ) is a Cauchy sequence. 

□ 

3.2.3 Existence and uniqueness of a solution of the mean-field equa- 
tions 

It is now easy to prove our main result, that the mean- field equations (|3ip or 
equivalcntly (|28|) arc wcll-poscd, i.e. have a unique solution. 

Theorem 3.10. For any nondegenerate kP -dimensional Gaussian random vari- 
able Zq, independent of the Brownian, and any initial process X such that 
X(to) = Zq, the map has a unique fixed point in (C({t , T], lR fcp )) to- 
wards which the sequence {^ n) (X)}£Li of Gaussian processes converges in law. 

Proof. Since C([t , T], R fep ) (respectively C([t a ,T] 2 , R fePxfep )) is a Banach space 
for the uniform norm, the Cauchy sequence fx n (respectively C n ) of proposi- 
tion [379] converges to an element \i of C([to, T], M. kp ) (respectively an element 
C of C([t ,T] 2 1 M, kPxkp )). Therefore, according to theorem O the sequence 

{^™' ) (^)}^Lo °f Gaussian processes converges in law toward the Gaussian pro- 
cess Y with mean function /1 and covariance function C . This process is clearly 
a fixed point of Th- 

Hence we know that there there exists at least one fixed point for the map Tk- 
Assume there exist two distinct fixed points Y\ and Y% of Tk with mean functions 
Hi and covariance functions Cj, i = 1,2, with the same initial condition. Since 
for all n > 1 we have .Ft™ (li) = Yi, i = 1,2, the proof of proposition 13.91 shows 
that — ^2 1 1 00 (respectively ||C" — C^|| ) is upperbounded by the product 
of a positive number a n (respectively b n ) with — (respectively with 

\\Cx — C^Hoo)- Since lim„_ ) . 00 a n = lim„^oo b n = and fif = fi iy i = 1,2 
(respectively C™ = Cj, i = 1,2), this shows that \i\ = fi 2 and C\ — C2, hence 
the two Gaussian processes Y\ and Y 2 are indistinguishable. 

□ 

Conclusion 

We have proved that for any non degenerate Gaussian initial condition Zq there 
exists a unique solution of the mean-field equations. The proof of theorem 13. 101 
is constructive, and hence provides a way for computing the solution of the 
mean-field equations by iterating the map Th defined in 13 . 21 starting from any 
initial process X satisfying X(to) — Zq, for instance a Gaussian process such as 
an Ornstein-Uhlcnbcck process. We build upon these facts in section [5] 

Note that the existence and uniqueness is true whatever the initial time to 
and the final time T. 
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4 Existence and uniqueness of stationary solu- 
tions 



So far, we have investigated the existence and uniqueness of solutions of the 
mean-field equation for a given initial condition. We are now interested in 
investigating stationary solutions, which allow for some simplifications of the 
formalism. 

A stationary solution is a solution whose probability distribution does not 
change under the flow of the equation. These solutions have been already in- 
vestigated by several authors (see 0, [34|). We propose a new framework to 
study and simulate these processes. Indeed we show in this section that under 
a certain contraction condition there exists a unique solution to the stationary 
problem. As in the previous section our proof is constructive and provides a 
way to simulate the solutions. 

Remark. The long-time mean-field description of a network is still a great 
endeavor in mathematics and statistical physics. In this section we formally 
take the mean-field equation we obtained and let to — > — oo. This way we obtain 
an equation which is the limit of the mean-field equation when to —* — oo. It 
means that we consider first the limit N — > oo and then to — ► — oo. These two 
limits do not necessarily commute and there are known examples, for instance 
in spin glasses, where they do not. 

It is clear that in order to get stationary solutions, the stochastic system has 
to be autonomous. More precisely, we modify assumptions 13 . 1 1 as follows 

Assumption 4.1. 



1. The matrixes L(t) and F(t), the input currents I(t) do not depend upon 
t. 

2. The real parts of the eigenvalues of L are negative: 

Re(A) < -X L A L > (48) 

for all eigenvalues A of L. 

3. The matrix F has full rank. 

Under assumption a of 14.11 the resolvent ^(i, s) is equal to e L ' t_s - ) . Under 
assumption 2 we only consider first-order system since otherwise the matrix L 
has eigenvalues equal to 0. We now prove the following proposition. 

Proposition 4.1. Under the previous assumptions we have: 

1. 



lim e L (*-*°) = 0, 

£q— >-— OO 
t 



S 



-(*-) || ds = f^ 



3 L T (t-s) 



ds = So 



oo 1 1 t 1 1 , def , 

"e Lu du = M L < oo, 

du = M l t < oo, 

oo 
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2. the process Y"/ = f t e L (* S ^F ■ dW s is well-defined, Gaussian and sta- 
tionary when to — ► — oo. 

Proof. The first point property follows from the fact that Re(A) < — Xl for all 
eigenvalues A of L. This assumption also implies that there exists a norm on 
R p such that 

|| e Lt|| < e -\ L t y t > Q) 

and hence 

lle"!^ < ke~ XLt Vt > 0, (49) 

for some positive constant k. This implies the remaining two properties. 

We now address the second point of the property. The stochastic integral 
Y t to = j ta e^-^F-dWs is well-defined Vt < T and is Gaussian with zero-mean. 
Its covariance matrix reads: 

rtht' 

s n to n,° = / e L(Mffr e LT ( t '-)d,. 



Let us assume for instance that t' < t and perform the change of variable 
u = t — s to obtain 

Under the previous assumptions this matrix integral is defined when to ~~ * 
— oo (dominated convergence theorem) and we have 

lim £ W ^ £^~>7~ = ( [ +C ° e^FF T e^ u du) e^^'^K (50) 

t ^-oo \Jt-t' J 

which is a well defined function of t' —t. □ 
The second point of propostion 14.11 guarantees the existence of process 

ri 



X (t) = / e L(f - !) F-<iW s . 

J — oo 



as the limit of the processes F t *° when to — * — oo. This process is a stationary 
distribution of the equation: 

<£X (f) =L-X (t)dt + F ■ dW t , (51) 

it is Gaussian, of mean E[Xo(t)] = and of covariance matrix S° is equal to 
T, Y * Yt defined by equation ([50]) and which is independent of t. 

We call long term mean-field equation (LTMFE) the implicit equation: 

V(f)= f e L (*- s )(u7-l + l)ds + X (t) (52) 

J — OO 
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where X is the stationary process defined by equation (|5Tj) and where U v (<) 
is the effective interaction process introduced previously. 

We next define the long term function JF stat : Mf(C((— oo, T], R p ) — ► 
.M+(C((-oo,T],R p ): 

e L( *- s) (uf -l + l) ds + X (t). 

-oo 



Proposition 4.2. The function J^at is well defined on Mi(C((— oo, T], R p ). 

Proof. We have already seen that the process Xo is well defined. The term 
f e L (*~ s )l<is = ( f* e L(t-s) ^ s ] J j s a j so we ii defined because of the as- 

J — oo w — oo y 

sumptions on L. 

Let X be a given process in A / J^"(C((— oo, T], R p ). To prove the proposition 
wc just have to ensure that the Gaussian process f_ e L (* _s )U^ ■ 1 ds is well 
defined. This results from the contraction assumption on L and the fact that 
the functions Sp are bounded. We decompose this process into a "long memory" 
term f_ e L ^ _s ^U^ • Ids and the interaction term from time t — 0, namely 

J 'e L ^ s 'Uf • Ids. This latter term is clearly well defined. We show that the 
memory term is also well defined as a Gaussian random variable. 

We write this term e Lt f e _Ls U? : • 1 ds and consider the second factor. 

J — oo s 

This random variable is Gaussian, its mean reads e Ls ^ u - s • 1 ds where 

p 



a=l...P 



A* s = I 2^J a ^[Sp{X {-s))]+I a 

The integral defining the mean is well-defined because of (|49|) and the fact that 
the functions Sp are bounded. A similar reasoning shows that the corresponding 
covariance matrix is well-defined. Hence the Gaussian process Jl^e-^Uf -Ids 
is well defined, and hence for any process X e A^^(C((— oo, T], R p ), the process 
^stat(X) is well defined. □ 

We can now prove the following proposition. 

Proposition 4.3. The mean vectors and the covariance matrices of the pro- 
cesses in the image of ^tat are bounded. 



Proof. Indeed, since E [-X" (t)] = 0, we have: 



e^-^u^ds 



<M L (H+\\I\U=HLT. 
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In a similar fashion the covariance matrices of the processes in the image of 
^stat are bounded. Indeed we have: 

E [^ stat (X) t ^ stat (X)f] = S°+ 
I I gLft-Odiag ( cr 2 a0 E [S (X ( Sl ))S (X (s 2 ))] I e^^) d Sl ds 2 , 

J — OO J — OO \ a J 

resulting in 

||E [^ tat (X)^ 8tat (X)F]L < ||S°|L +k 2 (^) 2 = Slt- 

□ 

Lemma 4.4. The set of stationary processes is invariant by J^tat- 

Proof. Since the processes in the image of J^tat are Gaussian processes, one 
just needs to check that the mean of the process is constant in time and that 
its covariance matrix C(s, t) only depends on t — s. 

Let Z be a stationary process and Y = jr stat (Z). We denote by fj, z the mean 
of the process Z a (t) and by C z (t — s) its covariance function. The mean of the 
process U Z p reads: 

mf /3 (*)=E[^(Z /3 (t))] = -== / ^(xje 2C ^ (0) dx 

^2^(0) -/B 

and hence does not depends on time. We note fi z the mean vector of the 
stationary process XJ Z ■ 1. 

Similarly, its covariance function reads: 

A z (t,s) = E[S (Z (t))S (Z (s))} = 

which is clearly a function, noted A z (t — s), of t — s. Hence XJ Z ■ 1 is stationary 
and we denote by C uZ (t — s) its covariance function. 
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It follows that the mean of Y t reads: 
= E[X (t)]+E 



(-: 



H*->) (i + uf ■ 1) ds 



f e Ht-s) (i + E[Uf -1]) ds 

J — OO 

e Lu dv) (I + f i Z ) 



Since we proved that E [Uf • l] = /j, z was not a function of s. 

Similarly, we compute the covariance function and check that it can be writ- 
ten as a function of (t — s). Indeed, it reads: 



C Y (t,s) 



[ [ e L(t - M) Cov(Uf • 1, Uf • l)e LT(s -' u) dudv + Cov{X (t), X {s)) 

J —oo J —oo 



rO 



C u (t - s + (u - v))e L v dudv + Cov(X (t),X {s)) 



oo J — oo 



since the process Xq is stationary. C Y (t, s) is clearly a function of t — s. Hence 
Y is a stationary process, and the proposition is proved. 

□ 

Theorem 4.5. The sequence of processes {J r s^at( x )}n'=o * s uniformly tight. 

Proof. The proof is essentially the same as the proof of theorem 13.41 since we 
can write 

T stllt (X)t=e u T stllt (X) + / e L <'->(U? -1 + I)ds + fe L « f -"»FdW s 

Jo Jo 

Fsta.t{X) t appears as the sum of the random variable -F s tat(^)o and the Gaussian 
process defined by f Q e L (*" s )(Uf • 1 + I) ds + f* e^'^ FdW s which is equal 



to T k {X) t defined in section for t = 0. Therefore F^t( x )t = -T^POt 
for t > 0. We have proved the uniform tightness of the sequence of processes 
{•^1 ( x )}'n-o m theorem 13.41 Hence, according to Kolmogorov's criterion for 
tightness, we just have to prove that the sequence of Gaussian random variables: 

^Stoo = ( f <S>l{-u){\3^ {X) ■ 1 + l)du + X (0) 



is uniformly tight. Since it is a sequence of Gaussian random variables, it is 
sufficient to prove that their means are bounded and their covariance matrices 
upperbounded to obtain that for any e > there exists a compact K e such that 
for any n € N, we have P(J r f t ™' ) t (X)o G K £ ) > 1 — e. This is a consequence of 
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proposition 14. 31 for the first random variable and of the definition of Xo for the 
second. By Kolmogorov's criterion the sequence of processes {J~^ t (X)} < ^' =0 is 
uniformly tight 

□ 

In order to apply theorem 13.11 we need to prove that the sequences of co- 
variance and mean functions are convergent. Unlike the case of to finite, this 
is not always true. Indeed, to ensure existence and uniqueness of solutions in 
the stationary case, the parameters of the system have to satisfy a contraction 
condition, and proposition 13.91 extends as follows. 

Proposition 4.6. If Al defined in (|4"5| satisfies the conditions ([53")) defined in 
the proof, depending upon kc (defined in l|45p). fco, ult and £lt (defined in 
proposition 14. 3p thcn the sequences of covariance matrix functions C n (t, s) and 
of mean functions /i n (t), s, t in [ioj^l & r e Cauchy sequences for the uniform 
norms. 

Proof. The proof follows that of proposition 13.91 with a few modifications that 
we indicate. In establishing the equation corresponding to (j44|) we use the fact 
that ||3>l(£, ii)!!^ < fce~ AL ( t_u ) for some positive constant k and all u, t, u < t. 
We therefore have: 



k 2 e -X L (t+s) 



D \ L (u+v) 



— oo J — oo 



Cov (U*« , Uf «)-Cov (U*- 1 , U*~ 



du dv 



The rest of the proof proceeds the same way as in proposition 13.91 Equations 
(1351) and H3 become: 



C n ^(t, s )-C n (t,s)\< 



Ke -\L{t+s) 



D \ L (u+v) 



[-oo.tvsp y/f{u,v) 

e \ L (u+v) 

[-oo,tvs] 2 y/f(u,v) 



C n (u,v)-O n - 1 (u,v)\\ dudv+ 



C n (u,u) - C" _1 (u,u)|| dudv+ 



-oo, tVs] 2 



e \ L (u+v) \\c n { u ,v) - C n - X {u,v)\\ dudv+ 



f e \L(u+v) \\ C n (u,u) - C^i^u)]^ dudv+ 

J[-oo,tVs] 2 



,\ L (u+v) \\..nr \ _ „n-l 



[-oo, tVs] 2 



M n (w) -/x" (u) I dudv), 
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and 

||/i" +1 (i)-/Li n (i)|| oo <Ke~ XL( - t+s H f e^f^ ||C n («,«)-C n - 1 (ti J u)|| oo dudv+ 

V J [—00, tVs] 2 

/ e A ^+")|| M >)-^-V)L^Y 

ii-octvs] 2 00 y 

for some positive constant K, function of k, kc (defined in (|45p ). and fco- 

Proceeding recursively until we reach C° and p° we obtain an upperbound 
for \\C n+1 (t,s) - C n (t,s) || x (respectively for ||^™ +1 (i) - ^"(i)^) which is the 
sum of less than 5™ terms each one being the product of K n , times 2/xlt or 
2Y,lt, times a 2n-dimensional integral /„ given by: 

Pi("i,wi) / ■• / Pn-i{u n -i,v n -i) 

[— oo.iVs] 2 \ J [— 00, ui V^i] 2 V J [— oo,u Tl „2V , u T1 -2] 2 

/ e Xh ^ Un+Vn ^ p n {u n ,v n )du n dv n \du n ^idv n -i\ ■ ■ ■ \duidvi, 

\ J[-oo,u„_ 1 V«„_ 1 ] 2 / J J 

where the functions pi (it,, Vi) , i = 1, • ■ ■ ,n are either equal to 1 or to 1 / \J 6(ui, Vi). 

It can be shown by straightforward calculation that each sub-integral con- 
tributes at most either 

A| lf Pi = 1 ° r 
in the other case. Hence we obtain factors of the type 

— 1 \ n ~ P / Hp, n -P / 1 \ (3n+p)/2 



MatJ [fr-tf) teJ *?• 

where < p < n. If A L < 1, (A L )( 3n +?>)/ 2 > Af™ and else (A l )( 3 "+p)/ 2 > xf /2 . 
Since (a/J) P < (\/1D we °btain the two conditions 



1 > A| > B^KKq or |a^ /2 > 5^KK and A L > l| (53) 



□ 



Putting all these results together we obtain the following theorem of exis- 
tence and uniqueness of solutions for the long term mean-field equations: 

Theorem 4.7. Under the contraction conditions (|53j) . the function T s t a t has a 
unique solution in A4^(C((— 00, T], R p ) which is stationary, and for any process 

X , the sequence {J r f™^ ( (X)}^L of Gaussian processes converges in law toward 
the unique fixed point of the function !F s tat- 
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Proof. The proof is essentially similar to the one of theorem 13.101 Indeed, the 
mean and the covariance matrixes converge since they are Cauchy sequences in 
the complete space of continuous functions equipped with the uniform norm. 
Using thcorcm l3.ll wc obtain that the sequence converges to a process Y which 
is necessarily a fixed point of J^stat- Hence we have existence of a fixed point for 
^-"stat- The uniqueness comes from the results obtained in the proof of propo- 
sition 14.61 The limiting process is necessarily stationary. Indeed, let X be a 
stationary process. Then for any n G N, the process .F st ™t(X) will be stationary 
by the virtue of lemma 14. A\ and hence so will be the limiting process which is 
the only fixed point of T s t a t . 

□ 

Hence in the stationary case, the existence and uniqueness of a solution 
is not always ensured. For instance if the leaks arc too small (i.e. when the 
time constants of the decay of the membrane potentials are too long) then the 
sequence can diverge or have multiple fixed points. 

5 Numerical experiments 
5.1 Simulation algorithm 

Beyond the mathematical results, the framework that we introduced in the 
previous sections gives us a strategy to compute numerically the solutions of 
the dynamic mean-field equations. Indeed, we proved in section [3] that under 
very moderate assumptions on the covariance matrix of the noise, the iterations 
of the map Th starting from any initial condition converge to the solution of 
the mean- field equations. 

This convergence result gives us a direct way to compute numerically the 
solution of the mean-field equations. Since we are dealing with Gaussian pro- 
cesses, determining the law of the iterates of the map amounts to computing 
its mean and covariance functions. In this section we describe our numerical 
algorithm in the case of the Simple Model of section 12.2.21 

5.1.1 Computing Tk- 

Let X be a P-dimensional Gaussian process of mean [i x = (fi x (t)) a =i...p and 
covariance C x = (C x Js, t)) a ,pe{i...p} ■ We fix a time interval {to = 0,T] and 
denote by Y the image of the process X under T\. In the case of the simple 
model, the covariance of Y is diagonal. Hence in this case the expressions we 
obtain in section [3] simply read: 

^(0)e-</^ + J* e-<*->/ T «(£j=i J af3 E [S (X (s))\ + I a (s))ds 
^(0)e-*/r«, + £e-M/T»I a { a )ds 
+ ELi J a p So e - (t " s)/T ° C Sp (x^f(s) + 4 (a)) Dxds. 



/£(*) = 
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where we denoted v% (s) the standard deviation of X a at time s, instead of 
C^ a (s,s). Thus, knowing v^(s),s £ [0,t] we can compute (J%(t) using a stan- 
dard discretization scheme of the integral, with a small time step compared with 
r a and the characteristic time of variation of the input current I a . Alternatively, 
we can use the fact that /i^ satisfies the differential equation: 



dt 



Ma 



P 

E 

/3=1 



+ 00 



= -— + Z^J a p / SplxJvf(t)+i$(t))Dx + I a (t), 



and compute the solution using a Runge-Kutta algorithm (which is faster and 
more accurate). Note that, when all the standard deviations of the process X are 
null for all time t £ [0, T] , we obtain a standard dynamical system. Nevertheless, 
in the general case, Vp(t) > for some /3's, and the dynamical evolution of (j% 
depends on the Gaussian fluctuations of the field X. These fluctuations must 
be computed via the complete equation of the covariance diagonal coefficient 
C^ a (t, s), which reads: 



CL(t, s) = e-(*+')/- [t£(0) + Zf^ ( e £ - l) 



/j=i 



Jo 



z {u+v)/T "A^Ju,v)dudv 



where: 



yjvf{u)vf{v)-Cf^u,vY Gf {u,v) 



X Sp (yy/v£(v) + $ («)) Dx Dy. 



Unless if we assume the stationarity of the process (see e.g. section I5T2|) . 
this equation cannot be written as an ordinary differential equation. We clearly 
observe here the non-Markovian nature of the problem: C* a (t,s) depends on 
the whole past of the process until time t V s. 



This covariance can be split into the sum of two terms: the external noise 



contribution C°£(t,s) = e -(*+ s )/ r ° 



(0) 



(eft-!) 



and the interac- 



tion between the neurons. The external noise contribution is a simple function 
and can be computed directly. To compute the interactions contribution to the 
standard deviation we have to compute the symmetric two-variables function: 



3 -(t+a)/ 7 



JO 



e^ +v)/Ta A^(u,v)dudv, 



from which one obtains the standard deviation using the formula 
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CZJt, ») = C°?(t, s) + J2 vlpH^it, s). 

To compute the function H x g (t,s), we start from t = and s = 0, where 
H x J0, 0) = 0. We only compute H x 3 (t, s) for t > s because of the symmetry. 
It is straightforward to see that: 



H x Jt + dt,s)=H x Jt,s) 



Ta 



+ D a c /3 (t,s)dt + o(dt), 



with 



Jo 



Hence computing H X p(t+dt, s) knowing H X g(t, s) amounts to computing D a p(t, s). 
Fix t > 0. We have D a p(t, 0) = and 



D x a (t, s + ds) = D^{t, s)(l - — ) + A^(t, s)ds + o{ds) 



ds. 



x 



This algorithm enables us to compute H x g (t, s) for t > s. We deduce H x g (t, s) 
for £ < s using the symmetry of this function. Finally, to get the values of 
H X g(t, s) for t — s, we use the symmetry property of this function and get: 



H^(t + dt,t + dt) = H x 3 (t, t) 



2dt 



+ 2D*p(t,t)dt + o(dt). 



These numerical schemes provide an efficient way for computing the mean 
and the covariance functions of the Gaussian process T\{X) (hence its proba- 
bility distribution) knowing the law of the Gaussian process X. The algorithm 
used to compute the solution of the mean-field equations for the general models 
GM1 and GMfc is a straightforward generalization. 



5.1.2 Analysis of the algorithm 

Convergence rate As proved in theorem 13.101 given Z a nondegenerate 
fcP-dimensional Gaussian random variable and X a Gaussian process such that 
X(0) = Z , the sequences of means and covariance functions computed theo- 
retically converge uniformly towards those of the unique fixed point of the map 
Tk- It is clear that our algorithm converges uniformly towards the real function 
it emulates. Hence for a finite N, the algorithm will converge uniformly towards 
the mean and covariance matrix of the process !F X {X). 

Denote by X f the fixed point of T k in Mt (C{[t , T], R fep )), of mean fi Xf {t) 

and covariance matrix C x - f (t,s), and by J- X (X) the numerical approximation 
of !F^(X) computed using the algorithm previously described, whose mean 

is noted /r 77 * ( x \t) and whose covariance matrix is noted C^fe W(t, a). The 
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uniform error between the simulated mean after N iterations with a time step 
dt and the fixed point's mean and covariance is the sum of the numerical error 
of the algorithm and the distance between the simulated process and the fixed 
point, is controlled by: 

||^ W - n x f 1^ + ||cW) - C x f IU = 0( {N + T) dt + i? Ar (fc max ) ) (54) 

where fc max = max(/c, k) and k and k) are the constants that appear in the 
proof of proposition 13.91 for the mean and covariance functions, and Rn(x) is 
the exponential remainder, i.e. Rn(x) = X^^at 2 '"/ 71 -'- 
Indeed, we have: 

|,/?P0 _ ^ f]]oo < ]]f f* ( x) ^ ( X) U + y ^ { x) _ ^ /||oo (55) 

The discretization algorithm used converges in O(dt). Let us denote by C\ 
the convergence constant, which depends on the sharpness of the function we 
approximate, which can be uniformly controlled over the iterations. Iterating 
the numerical algorithm has the effect of propagating the errors. Using these 
simple remarks we can bound the first term of the righthand side of (|55|) . i.e. 
the approximation error at the iVth iteration: 

||^ (X) -/^ W IU <C!iVdt 

Because the sequence of means is a Cauchy sequence, we can also bound the 
second term of the righthand side of (|55|) : 

oo 

H/^w < E \\^ +1[x) - ^ {x) u 

n—N 
oo , „ 

for some positive constant k introduced in the proof of proposition 13.91 The 
remainders sequence (R n {k)) n >o converges fast towards (an estimation of its 
convergence can be obtained using the fact that limsup fe ^ 00 (l/fc!) 1/ ' fe = by 
Stirling's formula) . 
Hence we have: 

||M^ (X) -A* X/ l|oo <C!Ndt + R N (k) (56) 
For the covariance, the principle of the approximation is exactly the same: 

_ cX,^ < llc ^ m C F»(X) lloo + ^(X) C X f | U 

The second term of the righthand side can be controlled using the same eval- 
uation by i?Ar(fc) where k is the constant introduced in the proof of proposition 
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13 . 01 and the first term is controlled by the rate of convergence of the approx- 
imation of the double integral, which is bounded by C 2 (TV + T) dt where C 2 
depends on the parameters of the system and the discretization algorithm used. 
Hence we have: 

||C^ Np0 - C Xf Woo < C 2 (TV + T - t ) dt + R N {~k) (57) 

The expressions ((56)) and (|5T[) are the sum of two terms, one of which is 
increasing with TV and T and decreasing with dt and the other one decreasing 
in TV. If we want to obtain an estimation with an error bounded by some e > 0, 
we can for instance fix TV such that max(R^(k), i?jv(&)) < e/2 and then fix the 
time step dt smaller than min( e/(2CiTV), e/(2C 2 (N + T- to)) )■ 



Complexity The complexity of the algorithm depends on the complexity of 
the computations of the integrals. The algorithm described hence has the com- 
plexity 0(TV(£) 2 ). 



5.2 The importance of the covariance: Simple Model, one 
population. 

As a first example and a benchmark for our numerical scheme we revisit the 
work of Sompolinsky and coworkers (34|. These authors studied the case of 
the simple model with one population (P = 1), with the centered sigmoidal 
function S(x) = tanh(gx), centered connectivity weights J = of standard 
deviation a = 1 and no input (J = 0,A = 0). Note therefore that there is 
no "noise" in the system, which therefore does not match the non degeneracy 
conditions of proposition ^ . 71 and of thcorcm l3.10l . This issue is discussed below. 
In this case, the mean equals for all t. Nevertheless, the Gaussian process is 
non trivial as revealed by the study of the covariance C(t, s). 



5.2.1 Stationary solutions 

Assuming that the solution of the mean-field equation is a stationary solution 
with C(t,s) = C(t — s) — C(t), Sompolinsky and his collaborators found that 
the covariance obeyed a second order differential equation : 

d2 ° dV « (58) 
~d^-~W (58) 

This form corresponds to the motion of a particle in a potential well and it is easy 
to draw the phase portrait of the corresponding dynamical system. However, 
there is a difficulty. The potential V q depends on a parameter q which is in 
fact precisely the covariance at r = (q = C(0)). In the stationary case, this 
covariance depends on the whole solution, and hence cannot be really considered 
as a parameter of the system. This is one of the main difficulties in this approach: 
mean-field equations in the stationary regime are self-consistent. 
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Nevertheless, the study of the shape of V q , considering q as a free parameter 
gives us some informations. Indeed, V q has the following Taylor expansion (V q 
is even because S is odd): 

V g (C) = ~C? + Ic' + oic 6 ) 

where A = (1 - g 2 J 2 (S') q ) and 7 = g J 2 g 6 (S^) 2 ), (4>) q being the average value 
of <j) under the Gaussian distribution with mean zero and variance q = C(0). 

If A > 0, i.e. when g 2 J 2 (S') 2 < 1, then the dynamical system (|55|) has a 
unique solution C(t) = 0,Vi > 0. This corresponds to a stable fixed point (i.e. 
a deterministic trajectory, = with no fluctuations) for the neural network 
dynamics. On the other hand, if g 2 J 2 (S'} 2 > 1 there is a homoclinic trajectory 
in (f58|) connecting the point q = C* > where V q vanishes to the point C = 0. 
This solution is interpreted by the authors as a chaotic solution in the neural 
network. A stability analysis shows that this is the only stabk0 stationary so- 
lution (3i]. 

The equation for the homoclinic solution is easily found using energy con- 
servation and the fact that V q (q) — and ^f(q) = 0. One finds: 

At the fourth order in the Taylor expansion of V q this gives 




Though A depends on q it can be used as a free parameter for interpolating 
the curve of C(r) obtained from numerical data. 



5.2.2 Numerical experiments 

This case is a good benchmark for our numerical procedure since we know 
analytically the solutions we are searching for. We expect to find two regimes. 
In one case the correlation function is identically zero in the stationary regime, 
for sufficiently small g values or for a sufficiently small q (trivial case). The 
other case corresponds to a regime where C(t) > and C(r) — > has t — > +00 
("chaotic" case). This regime requires that g be sufficiently large and that q 
be large too. We took t q = 0.25, a aa = 1. For these values, the change in 
dynamics predicted by Sompolinsky and collaborators is g c = 4. 

In sections [3] and [4] we have introduced the assumption of non-degeneracy 
of the noise, in order to ensure that the mean-field process was non degenerate. 
However, in the present example, there is no external noise in the evolution, 
so we can observe the effects of relaxing this hypothesis in a situation where 

7 More precisely, this is the only minimum for the large deviation functional. 
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Figure 3: Numerical solution of the mean-field equation after 14 iterations in the 
chaotic case (g = 5). We clearly see the numerical instabilities in the no- noise 
case, which do not exist in the low-noise case. 

the results of proposition 13.71 and of theorem 13.101 cannot be applied. First, we 
observed numerically that, without external noise, the process could become 
degenerate (namely some eigenvalues of the covariance matrix C a (t, s) become 
very small and even vanish.). This has also an incidence on the convergence of 
the method which presents numerical instabilities, though the iterations leads 
to a curve which is well fitted by the theoretical results of Sompolinsky et al. 
(see Fig. [3]) . The instability essentially disappears if one adds a small noise. 
But, note that in this case, the solution does not match with Sompolinsky et 
al. theoretical calculation (see Fig. [3]). 

Modulo this remark, we have first considered the trivial case corresponding 
to small g values. We took g = 0.5 and T = 5. We choose as initial process the 
stationary Ornstein-Uhlenbeck process corresponding to the uncoupled system 
with A = 0.1. We drew fi a (0) randomly from the uniform distribution in [—1, 1] 
and v a (0) randomly from the uniform distribution in [0, 1]. 

Starting from this initial stationary process, we iterated the function T\. 
Then, during the iterations, we set s a = in order to match the conditions im- 
posed by Sompolinsky and coworkers. We observe that the method converges 
towards the expected solution: the mean function converges to zero, while the 
variance v(t) decreases exponentially fast in time towards a constant value cor- 
responding to the stationary regime. This asymptotic value decreases between 
two consecutive iterations, which is consistent with the theoretical expectation 
that v(t) = in the stationary regime of the trivial case. Finally, we observe 
that the covariance C(t — s, s) stabilizes to a curve that does not depend on s 
and the stationary value (large t — s) converges to zero. 

We applied the same procedure for g = 5 corresponding to the "chaotic" 
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regime. The behavior was the same for but was quite different for the 
covariance function C(t, s). Indeed, while in the first case the stationary value of 
v(t) tends to zero with the number of iterations, in the chaotic case it stabilizes 
to a finite value. In the same way, the covariance C(t — s, s) stabilizes to a 
curve that does not depend on s. The shape of this curve can be extrapolated 
thanks to Sompolinsky et al. results. We observe a very good agreement with 
the theoretical predictions with a fit fn(x) — cos h(b°3:-g)) ' corresponding to the 
fourth expansion of V q . Using a 6-th order expansion of V q (x) = ^x 2 + jx 4 + ^x 2 
gives a fit f 6 (x) = coshl £ x _ s)) u , = 1 ^ — , where p,K,X are explicit 

+ co B h2(A(x-6)) 

functions of a, 6, c, we obtain a slightly better approximation. 

5.3 Mean-field equations for two populations with a neg- 
ative feedback loop. 

Let us now present a case where the fluctuations of the Gaussian field act on the 
dynamics of (i a [t) m a non trivial way, with a behavior strongly departing from 
the naive mean-field picture. We consider two interacting populations where the 
connectivity weights are Gaussian random variables J a p = Af(J a p, <J a p = 1) for 
(a, 0) € {1, 2} 2 . We set Sp(x) = tanh(.g.x) and I a = 0, s a = 0, a = 1, 2. 

5.3.1 Theoretical framework. 

The dynamic mean-field equation for /i a (t) is given, in differential form, by: 

Let us denote by G a (fi, v(t)) the function in the righthand side of the equal- 
ity. Since S is odd, S(y/vp(t)x)Dx — 0. Therefore, we have G a (0, v(t)) = 
whatever v(t), and hence the point [i\ = 0, fi2 — is always a fixed point of this 
equation. 

Let us study the stability of this fixed point. To this purpose, we compute 
the partial derivatives of G a (fi, v(t)) with respect to up for (a,/3) G {1,2} 2 . We 
have: 



^(ji,v(t)) = -^-+9J«0 J (^l-tanh 2 \ J^x + fip(t)jJ Dx, 
and hence at the point fj,\ = 0, [12 = 0, these derivatives read: 

^(0,«(t)) = + gj a ph(vp(t)), 

where h(vp(t)) = 1 — tanh 2 (yjvp (t)x)Dx. 
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In the case v a (0) = 0, J = 0, s a = 0, implying v a (t) = 0, t > 0, the equation 
for /x a reduces to: 



/3=1 



eft 



which is the standard Amari-Cohen-Grossberg-Hopfield system. This corre- 
sponds to the naive mean-field approach where Gaussian fluctuations are ne- 
glected. In this case the stability of the fixed point fi = is given by the sign 
of the largest eigenvalue of the Jacobian matrix of the system that reads: 



h \ , ( Ju Ji 



tl 



-i J21 J. 



12 



'22 



For the sake of simplicity we assume that the two time constants r a are equal 
and we denote this value r. The eigenvalues are in this case — i + gX, where A 
are the eigenvalues of J and have the form: 



, Ju + J22 ± \/{Jii - J22) 2 +4J12J21 
Ai, 2 2 ■ 

Hence, they are complex whenever Jyi J 21 < — (Jli — J22) 2 /4, corresponding 
to a negative feedback loop between population 1 and 2. Moreover, they have 
a real part only if Ju + J22 is non zero (self interaction). 

This opens up the possibility to have an instability of the fixed point (fx = 0) 
leading to a regime where the average value of the membrane potential oscillates. 
This occurs if Ju + J22 > and if g is larger than: 



9c 



t(Jh + J 22 ) 



The corresponding bifurcation is a Hopf bifurcation. 

The situation is different if one takes into account the fluctuations of the 
Gaussian field. Indeed, in this case the stability of the fixed point \i = de- 
pends on v(t). More precisely, the real and imaginary part of the eigenvalues 
of DG(Q,v(t)) depend on v(t). Therefore, the variations of v(t) act on the 
stability and oscillations period of v(t). Though the evolution of fi(t),v(t) are 
coupled we cannot consider this evolution as a coupled dynamical system, since 
v(t) = C(t, t) is determined by the mean-field equation for C(t, s) which cannot 
be written as an ordinary differential equation. Note that we cannot assume 
stationarity here, as in the previous case, since n(t) depends on time for suf- 
ficiently large g. This opens up the possibility of having complex dynamical 
regimes when g is large. 

5.3.2 Numerical experiments 

We have considered the case Ju = J22 = 5,t = 0.1 giving a Hopf bifurcation 
for g c = 2 when J = (fig. [4|. The trajectory of ni(t) and vi(t) is represented 
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in Figure [Hin the case g = 3. When J = 0, fii(t) presents regular oscillations 
(with non linear effects since g = 3 is larger than the critical value for the 
Hopf bifurcation, g c = 2). In this case, the solution v\(t) = is stable as seen 
on the figure. When J ^ the Gaussian field has (small) fluctuations which 
nevertheless strongly interact with the dynamics of n\ (t) , leading to a regime 
where fii(t) and vi(t) oscillate periodically 



1 1 1 


mu(t), J=0 - 
mu(t), J=2, n-100 

v(t), J=0 □ 
v(t) x 100, J=2, n=100 




■ 






RyJ\ AA A AA\ 


ll/A \AJf IKS U 



2 4 6 8 10 12 14 16 18 20 



Figure 4: Evolution of the mean fii(t) and variance vi(t) for the mean-field of 
population 1, for J = and J = 2, over a time window [0, 20] . n is the number of 
iterations of T\ defined in section [3] This corresponds to a number of iterations 
for which the method has essentially converged (up to some precision). Note 
that t>i(i) has been magnified by a factor of 100. Though Gaussian fluctuations 
are small, they have a strong influence on jJL\(t). 



6 Discussion 

The problem of bridging scales is overwhelming in general when studying com- 
plex systems and in particular in neuroscience. After many others we looked 
at this difficult problem from the theoretical and numerical viewpoints, hoping 
to get closer to its solution from relatively simple and physically/biologically 
plausible first principles and assumptions. One of our motivations is to better 
understand such phcnomcnological neural mass models as that of Jansen and 



Rit [27 1 



We consider several populations of neurons and start from a microscopic, 
i.e. individual, description of the dynamics of the membrane potential of each 
neuron that contains four terms. 

The first one controls the intrinsic dynamics of the neuron. It is linear in 
this article but this assumption is not essential and could probably be safely 
removed if needed. 
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The second term is a stochastic input current, correlated or uncorrelated. 
The corresponding noise model is very rich, depending on the degree k of 
smoothness of the g-shapes. It features integrated Brownian noise up to or- 
der k — 1. 

The third term is a deterministic input current, and the fourth one describes 
the interaction between the neurons through random connectivity coefficients 
that weigh the contributions of other neurons through a set of functions that are 
applied to their membranes potentials. The only hypothesis on these functions 
is that they are smooth and bounded, as well as their first order derivative. The 
obvious choice of sigmoids is motivated by standard rate models ideas. Another 
appealing choice is a smooth approximation to a Dirac delta function thereby 
opening a window on the world of spiking neurons. Thus, the model presented in 
this paper is more general than the instantaneous rate model that is underlying 
Ermentrout's voltage-based model [3] even though we have not explored this 
avenue. 

We then derive the mean-field equations and provide a constructive and new 
proof, under some mild assumptions, of the existence and uniqueness of a so- 
lution of these equations over finite and infinite time intervals. The key idea 
is to look at this mean-field description as a global problem on the probability 
distribution of the membranes potentials, unlike previous studies. Our proof 
provides an efficient way of computing this solution and our numerical exper- 
iments show a good agreement with previous studies. It is interesting to note 
that a sufficient condition for the convergence of our algorithm is related to the 
previously mentioned noise model. We prove that if the noise matrix F is full 
rank, with bounded eigenvalues, then the algorithm is in general convergent. 
An important fact to note is that the solutions of the mean-field equations that 
we construct are fundamentally non-Markovian eliminating the need for such 
approximations as the introduction of the q parameter summarizing the whole 
history of the non-Markovian process, see below. 

In the case where the nonlinearities are chosen to be sigmoidal our results 
shed a new light on existing neural mass models. Indeed, as shown in section 
I2.2.H these appear as approximations of the mean-field equations where the 
intricate but fundamental coupling between the time variations of the mean 
membrane potentials and their fluctuations, as represented by the covariancc 
functions, is neglected. 



An alternative approach has been recently proposed by Chizhov and col- 
laborator i@ El. The approach of these authors consists in reducing the 
large number, N, of degrees of freedom of the neural assembly by construct- 
ing a probability density p on the phase space of neurons states in the limit 
N — > oo. This is a non rigorous approach where the evolution equations for 
p are heuristically derived. Especially, it is assumed that p depends on two 
parameters only: the current time t and the time elapsed since the last spike 
t* . Under these assumptions the initial phase space of neurons states is mapped 



'We thank one of the referees for pointing out these references to us. 
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to a two dimensional space t, t* , while p(t,t*)dt characterizes the fraction of 
neurons which have fired in the time interval [t — t*,t — t* + dt]. Therefore, 
this approach intrinsically holds for integrate and fire neurons models where 
the neuron's membrane potential history is summarized by the last spike time, 
when it is reset to a constant value. As noticed by these authors, this allows 
to circumvent the main problem in mean-field approaches for firing rate mod- 
els, that we also discuss in the present paper: When using mean-field theory 
to characterize stationary regimes, one needs to introduce ad hoc parameters 
(see e.g. the parameter q introduced in section 15.2. ip summarizing the whole 
history of the non-Markovian process. Introducing a "history cut-off" while 
reseting the membrane potential to a constant value indeed removes this diffi- 
culty Therefore, it might be interesting to compare our approach in the case of 
integrate- and- fire models (see above remark on the choice of the nonlincarity), 
to the approach of Chizov and collaborators. This could provide some rigorous 
basis for their analysis and allow to elucidate the role of field fluctuations which 
does not appear explicitely in the probability density approach. 

7 Conclusion and further work 

On more general grounds, our goal is now to extend the present work in several 
directions. 

Bifurcations analysis of the dynamic mean-field equations. From the 
present analysis, and as shown in the simple examples of section [5l the meso- 
scopic dynamics of the average membrane potential of a neurons population can 
be really different from the classical phcnomenological equations a la Jansen-Rit 
if one includes the non-Markovian fluctuations of the interaction fields, which 
summarize the cumulative effects of the nonlinear interactions of a given neuron 
with the bulk of other neurons. Jansen-Rit equations are commonly used in the 
neuroscience community either to anticipate the dynamics of local field potential 
in relation with imaging (Optical Imaging, MEG-EEG), or to understand neuro- 
physiological disorders such as epilepsy. Bifurcations analysis of these equations 
reveal dynamical regimes that can be related to experiments (23j . They can be 
generalised using more accurate neural models [38j . Is there any need to gener- 
alize these equations, that we claim to be incomplete, while people commonly 
use them with some satisfaction? Are there new phenomena, experimentally 
accessible, that can be exhibited by the generalised mean-field equations and 
that do not appear in the naive ones? These are obviously important questions 
that we intend to address in the near future. On mathematical grounds, the 
goal is to make a bifurcation analysis of the map T on the space of trajecto- 
ries, introduced in the present paper. Do any new salient dynamical regimes 
appear? If such regimes exist, the goal will be, on experimental grounds, to 
interact with experimentalists in order to see in which conditions such a regime 
can be exhibited, and what are its implications on cortical columns dynamics 
or function. 
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Investigations of non stationary regimes. As discussed in this paper, and 
as is well-known in the physicists' community (especially spin-glasses commu- 
nity) , the dynamic mean-field approach raises serious difficulties as far as one is 
trying to describe stationary dynamics. On technical grounds, this relies on the 
non-commutativity of the two limits N — > oo and t — > oo already discussed in 
(35l | . As a result, one is led to introduce ad- hoc phenomenological parameters, 
depending on initial conditions, that can be determined in statistical physics 
models where the distribution of equilibria is known (Gibbs distribution), using 
sophisticated techniques such as the replica "trick" [20|. For spin-glasses it is 
only in the high temperature regime that a simple solution to this problem is 
known. This restriction also appears in the present paper, where the existence 
and uniqueness of a stationary solution is proved only for low values of the gain 
parameter g (which plays a role similar to the inverse temperature). However, 
we are not so much interested in stationary dynamics, since brain processes are 
ultimately non stationary. Our approach, valid for any finite time T, opens up 
the possibility to characterize mean-field equations in transient regimes, with 
an analysis strategy that can moreover be easily implemented. To the best of 
our knowledge, this type of techniques has never been used in the statistical 
physics community, where iterations on space trajectories are not in the stan- 
dard toolbox. Therefore, our work could allow the (numerical) investigation of 
cortical columns submitted to non stationary inputs, with strong implications 
on ncuroscience. 



Extension to a larger class of models. A very challenging question is the 
application of this theory to spiking neurons models. We have briefly mentioned 
in section[6]that this may be possible through the use of non-sigmoidal functions 
in the interaction terms. This idea could be applied to the analysis of Integrate 
and Fire models with conductance based synapses, which constitute good models 
of spiking neurons. As discussed at the end of section the analysis of the 
mean-field equations could be simplified by the fact that memory is reset after 
a neuron fires. There is however a need to characterize parameter space reg ions 



where neurons can take an arbitrary large time to fire for the first time [lll.ll2|. 
This is the main obstruction in the application of our theory to this type of 
models. 
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A Identification of the mean-field equations 

Ben-Arous and Guionnct studied from a mathematical point of view the prob- 
lem of finding a mean-field description of large networks of spin glasses. They 
obtained using different methods of stochastic analysis a weak limit of the law 
of a given spin and proved their independence. 

Our equations do not directly fit in their study: indeed, the spin intrinsic 
dynamics is nonlinear while the interaction is linear, and everything in done in 
dimension one. Nevertheless, their proof extends to our case which is somehow 
more simple. For instance in the case of the Simple Model with one population, 
we can readily adapt their proof in our case. More precisely, let P — 1, the 
equation of the network reads: 

N 

rdVl = (-V t j + J2 J»S(Vl)) dt + adWl 

i=l 

In this case, we define for X £ Mi(C([to, T], R) the effective interaction 
term (t/ t x ) which is the effective interaction process defined in 12.21 i.e. the 

Gaussian process of mean J a pE[S(Xt)] and of covariance: Cov(lJ*, C/ x 

a 2 a(3 E[S(X t )S(X s )}. 

Let us note V the law of the membrane potential when there is no in- 
teraction (it is an Ornstein-Ulhcnbeck process), and the empirical measure 
V N = -h Yli=i ^v*- We can prove that under the probability distribution aver- 
aged over the connectivities, see below, the empirical measure satisfies a large 
deviation principle with good rate function H defined as in [24[ ■ Using this large 
deviation result, we can prove annealed and quenched tightness of the empirical 
measure, and finally its convergence towards the unique process where the good 
rate function H achieves its unique minimum, which is defined by the property 
of having a density with respect to V and whose density satisfies the implicit 
equation: 



8 


exp < 


I 









T 

Q\2 



(U?Ydt 



(59) 



where £ denotes the expectation over the effective interaction process 

We can also prove following the steps of Ben-Arous and Guionnet in 0] that 
there exists a unique solution to this equation, and that this solution satisfies 
the nonlinear nonmarkovian stochastic differential equation: 
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rdVt = -V t dt + dB t 



dB t = dW t + f* dB s £ U®U?- 




(60) 



Law of (V) = Q, law of (V ) = Z a 



which can also be written as our mean-field equation, averaged on the con- 
nectivities (see 0]). More precisely, let L v be the law of the solution of the 
equation: 



which is exactly equation (|33[) . They prove that V satisfies the nonlinear equa- 
tion: 



This result is probably extendable to the multi-population case using the mul- 
tidimensional Girsanov's theorem, but the corresponding mathematical devel- 
opments are out of the scope of this paper. 

B The resolvent 

In this appendix we introduce and give some useful properties of the resolvent 
$i of a homogeneous differential equation 



where L : [t ,T] -► M PxP (or (—00, T] -► M PxP ) is C°. 

Definition B.l. The resolvent of (j6"T]) is defined as the unique solution of the 
linear equation: 



rdVt = -V t dt + dW t + UYdt 
Law of Vq = Zq 



V = £{L V ) 



— = L(t)x(t) x(t ) = x & R p , 



(61) 




(62) 



where Idp is the P x P identity matrix. 



Proposition B.l. The resolvent satisfies the following properties: 



(i). $ £ (t + *, to) = $ £ (t + a, t). $x,(t, to) 



(ii). $i(t, to) is invertible of inverse 4>^(t ,t) which satisfies: 




(63) 
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(iii) . Let || || be a norm on Mp x p and assume that ||L(t)|| < on [to,T]. 

Then we have: 

||*i(*,*o)|| <e^l 4 - to l Vt€[t ,T] (64) 
Similarly, if ||L T (i)|| < k L r on [iojT 1 ] we have: 

\\^l(t,t )\\ < e k ^-^ Vte[t ,T} (65) 

(iv) . We have 

det$ L (Mo) = exp / TrL(s)ds 

J t 

Proof. The properties (i) and (ii) are directly linked with the property of group 
of the flow of a reversible ODE. (iii) is an application of Gronwald's lemma, (iv) 
is obtained by a first order Taylor series expansion. □ 

Theorem B.2 (Solution of an inhomogeneous linear SDE). The solution of the 
inhomogeneous linear Stochastic Differential Equation: 

idX t = (L(t)X(t)+I(t))dt + F(s)dW s 

\x to = x 

can be written using the resolvent: 

X t = $ L (t,t )X + f $ L (t,s)I(s)ds+ [ $ L (s,t)F(s)dW 3 (67) 

Proof. Pathwise (strong) uniqueness of solution directly comes from the results 
on the SDE with Lipschitz coefficients (see e.g. [HI, Theorem 2.5 of Chapter 
5]). It is clear that Xt„ = Xq. We use Ito's formula for the product of two 
stochastic processes to prove that the process (|67p is solution of equation (f6"6")) : 

dX t = (h(t)$ L (t, t Q )X Q + <f> L (t, t)I(t) + jf L(t)$i(t, s)I(s) da) dt 
+ $ L (t,t)F(t)dW t + f Ut)$ L (s,t)F(s)dW s dt 

J to 

= (h(t)h L (t,t )X + f <f> L (s,t)I(s)ds+ [ ® L (s,t)F(s)dW s ]+I(t))dt 

V L Jtn Jtn J ' 



+ F(t)dW t 

= (L(t)X(t) + 1(f)) dt + F(t)dW t 
Hence the theorem is proved. □ 
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C Matrix norms 



In this section we recall some definitions on matrix and vector norms. Let M n xn 
be the set ofnxn real matrices. It is a vector space of dimension n 2 and the 
usual L p norms 1 < p < oo can be defined. Given L £ A4 nxn , we note ||L||p 
the corresponding norm. Given a vector norm, noted || ||, on ]R n the induced 
norm, noted j| ||, on A^ n xn is defined as 

NT II H La; ll 

l|L|| = sup — — 
sen™, IMI<i IFII 

Since M nxn is finite dimensional all norms are equivalent. In this article we 
use the following norms 

W- ll L lloo =maxiE"=i \Lij\. 

(ii) - II L IIL = max M \ L ij\ 

(iii) . ||L|| 2 = sup^gj^,, ||x||,<i ^x\f 2 • This so-called spectral norm is equal to 

the square root of the largest singular value of L which is the largest 
eigenvalue of the positive matrix L T L. If L is positive definite this is its 
largest eigenvalue which is also called its spectral radius, noted p(L). 



D Important Constants 

Table [T] summarizes some notations which are introduced in the article and used 
in several places. 



Constant 


Expression 


Defined in 


A 4 


max Q J2p \J<*p\ Halloo 


lemma 13.21 
equation (IBUl) 


'2 

^max 


max Q Y,g a ap \\ s p\\'L 


lemma 13.21 




mm a< p al 


lemma 13.21 




e w-t ) m [Zq]11oq + (/i + /max)(T _ to)] 


lemma 13.51 


^max 




lemma 13.51 




\ ■ 


lemma 13.61 


K 


Arain-y/A^ A^ in (T - t ) 


proof of lemma 13.81 


kc 


maXaEja^/jll^lloo S 'f} 

oo 


proposition 13.91 
equation fi5l) 


\l 




equation (|48j) 



Table 1: Some important quantities defined in the article. 
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E Proof of lemma 13.51 

Lemma E.l. The following upppcrbounds are valid for all n > 1 and all 

a, t g [t ,T]. 

||/i n (*)IL < e^ (T ^ 0) [||E [^o]IL + (M + W)(T - *b)] = /W, 

s)^ < c (^+**T)Cr-to) [ p(E *,) + aLx(t _ ^ + CTLx(r _ io)2] def Sm 

where /i and <7 max are defined in lemma [3T21 A^ ax is defined in 13. f I 

Proof. The first inequality follows from taking the infinite norm of both sides 
of equation (|41[) and using assumption 1. in 13.11 and equation (|64[) . lemma [3~2l 
and assumption 3. in 13.11 

The second inequality follows from taking the infinite norm of both sides 
of equation (j4"2")) and using assumption 1. in 13.11 and equations and 
lemma 1321 and assumption 2. in 13.11 □ 



F Proof of lemma 13.61 

Lemma F.l. For all t <G [to, T] all a = 1, • • • , fcP, and n > 1, we have 

^ AminA^ju = fc > 0, 

where A m i n is the smallest singular value of the symmetric positive definite 
matrix to)&i,(t, t ) T for t £ [to,T] and A^;° is the smallest eigenvalue of 

the symmetric positive definite covariance matrix T, z ° . 

Proof. C£a(t,t) is larger than ($E,(i, t )Y, z °<fr L (t, t ) T ) aa which is larger than 
the smallest eigenvalue of the matrix $z,(t, to)T, z °^L(t, to) T . This smallest 
eigenvalue is equal to 

. x T $ L (t,i )£ Zo $L(Mo) T 2: 
mm = = 

||x||<l X 1 X 

min f x T <S> L (t, t )-Z z °$ L {t, t ) T x x T $ L (t, to)$ L (t, t ) T X 

IHI<i\ a; T $ L (t,t )$L(t,io) T a; x T x 

. x T $ L (t,t )Z z °<l> L (t,t Q ) T x . x T <S> L (t,t Q )<l> L (t,t ) T x 
mm = ; r ; r-= mm = . 

NI<1 x T $ L (t,t )<f> L (t,t ) T x \\x\\<i x T x 

In the last expression the first term is larger than the smallest eigenvalue 
A^;° of the matrix S z ° which is positive definite since we have assumed the 
Gaussian random variable Zq nondegencrate. The second term is equal to the 
smallest singular value A m i n of the matrix &L{t, to) which is also strictly positive 
since $z,(f, to) is invertible for all t € [to,? 1 ], see appendix [Bj 

□ 
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G Proof of lemma 13.7 

Lemma G.l. For all a = 1, ■ • ■ , kP and n > 1 the quantity C™ a (s, s)C™ a (t, t) — 
C™ Q (t, s) 2 is lowerbounded by the positive symmetric function: 

ft( s i\ — £ 1^ — «l X 2 \ T,Z ° \ r 

where is the strictly positive lower bound, introduced in l3.11 on the singular 
values of the matrix F(u) for u £ [to, T]. 

Proof. We use equation (|^|) which we rewrite as follows, using the group prop- 
erty of the resolvent <&l : 

C n+1 {t,s) = $ L (t,i )^£ z ° + <I> L (t Q ,u)F(u)F(u) T <Z> L (t ,u) T du+ 

jf jf $ L (i ,u)Cov(u*»,U*") $ L (t ,t.) T ^d^$ L ( s ,t ) T . 

We now assume s < t and introduce the following notations: 

A(s) = Z z °+f t S o <S> L (t ,u)F(u)F(u) T <S> L (t 0: u) T du 
B(s,t) = f*$L(to,u)F(u)F(u) T $ L (to,u) T du 
a(t,s) = /;j/ o $ i (t , u )Cov(u^,U^)$ L (t ,f) T du^ 

Let e a , a = 1, • • • , fcP, be the unit vector of the canonical basis whose coordi- 
nates are all equal to except the ath one which is equal to 1. We note E a (t) 
the vector to) T e a . We have, dropping the index n for simplicity: 

C aa (t,s) = E a (t) T (A(s) + a(t, s)) E a (s) 
Coals, a) = E a {s) T {A{s) +a(s,s))E a (s) 
C aa {t,t) = E a (t) T (A(s)+B(s,t)+a(t,t))E a {t). 

Note that the last expression does not depend on s, since A(s) + B(s, t) = A(t), 
which is consistent with the first equality. The reason why we introduce s in 
this expression is to simplify the following calculations. 

The expression C aa (s, s)C aa (t, t)—C aa (t, s) 2 is the sum of four sub-expressions: 

£i(a,t) = (E a ( s ) T A( s )E a (s)) (E a (t) T A(s)E a (t)) - (E a (t) T A(s)E a (s))\ 
which is greater than or equal to because A(s) is a covariance matrix, 

S 2 (s,t) = (E a (s) T a(s,s)E a (s)) (E Q (t) T a(t,t)E Q (t)) - (E a {t) J 'a(t, s)E a {s)) 2 , 

which is also greater than or equal to because a(t, s) is a covariance matrix 
function, 

S 3 {s,t) = (E a {s) T A(s)E a {s)) (E a {t) T a(t,t)E a {t)) + 

(E a (t) T A(s)E a (t)) (E a (s) T a(s,s)E a (s)) - 

2 (E a (t) T A(s)E a (s)) (E a (t) T a(t,s)E a (s)) 
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Because a(t, s) is a covariance matrix function we have 

E a (t) T a(t,t)E a (t) + E a (s) T a(s,s)E a (s) - 2E a (t) T a(t, s)E a (s) > 0, 
and , as seen above, £2(3, t) > 0. Because £i(s, t) > we also have 



- ^ E a (s) T A(s)E a (s)^ E a {t) T A(s)E a (t) < E a (tf A(s)E a (s) < 

^E a ( s ) T A( s )E a ( s )^E a (t) T A( s )E a (t), 

and, as it can be readily verified, this implies ^(s, t) > 0. 

Therefore we can lowerbound C aa (s, s)C aa (t, t) — C aa (t,s) 2 by the fourth 
subexpression: 

C aa {s,s)C aa (t,t)-C aa {t,s) 2 > (E a {s) T A{s)E a {s)) (E a (tf B(s,t)E a {t)) + 
(E a ( S ) T a{s,s)E a {s)) (E a {t) T B{s,t)E a (t)) > 

(E a (s) T A( S )E a (s)) (E a (t) T B(s,t)E a (t)) , 

since B(s, t) and a(s, s) are covariance matrixes. We next have 

p f.^i/,,^ E a (s) T A{s)E a (s) e^ L {s,t Q )^ L (s,to) T e a 
E a (s) A(s)E a (s)= Ea{s)TEa{s) ^ , 

by definition of E a (s). Therefore 

E a (sfA(s)E a (s) > A^A*L (s ' io) * l(Mo)T > A^A min , 

where A^ in is the smallest eigenvalue of the symmetric positive matrix C . Sim- 
ilarly we have 

E a (t) T B( S ,t)E a (t) > Afjf ] X min . 
Let us write T(u) = F(u)F(u) T . We have (assumptions 13 . 1 1 ) : 

x B(s,t) f* x T $ L (to,u)T(u)$ L (t ,u) T x 
A, rl in = nun / = du = 

iHi<iy a x T x 

min f* x T <$> L (t , u)T{u)^ L (t 0j u) T x x T $ L (t , u)$ L (t , u) T x > 
\\x\\<ij a x T $ L (t ,u)$ L (t ,u) T x x T x 

min / x T $ L (t , u)T(u)<f> L (t , u) T x x T <S> L (t ,u)<P L (t ,u) T x \ ^ > 
VM\<i\ x T $ L (t ,u)<i>L(t ,u) T x x T x ) U ~ 

(t — s)A m i n A r 



nun • 



Combining these results we have 

□ 
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H Proof of lemma 13.81 



Lemma H.l. The 2n-dimcnsional integral 



I n = / p±(ui,vi)l / " / Pn-i(u n -i,V n -i) 

J[t .tWs] 2 \ J[t ,«iVt;i] 2 \J[i ,ii„_2Vu„_2] 2 

/ p n (u n ,v n )du n dv n \du n - 1 dv n - 1 \---\du 1 dv 1 , 
\./[i ,u, l _ 1 vu„_ 1 ]2 I J J 



where the functions pi(u%,Vi), i = 1, • • • ,n are either equal to 1 or to 1 / \J 0(ui, m) 
(the function 8 is defined in lemma [377)1 . is upperbounded by k n /(n — 1)1 for 
some positive constant k. 

Proof. First note that the integral is well-defined because of lemma 13.71 Sec- 
ond, note that there exists a constant K such that K/ \/9(u, v) > 1 for all 

(u,v) € [i ,(Vsf, i.e. K = A min ^ A^f° A^ in (T — to)- Therefore the integral is 
upperbounded by Kfi , where Kq = max(l, K) times the integral obtained when 
Pi(ui, Vi) = 1/ \J\ui — Vi \ for alii = 1, • • • , n. Let us then consider this situation. 
Without loss of generality we assume to = 0. The cases n = 1, 2, 3 allow one to 
understand the process. 

f dudv .„„. 
A < #o / 7= = . (68) 

J[0,tVs] 2 v l u - v \ 

Let us rotate the axes by — j by performing the change of variables 

U + V 

V - u 

Using the symmetry of the integrand in s and t and the change of variable, 
the integral in the righthand side of flSSJ) is equal to (see figure [5]): 

2 l/i J Ju VU Jo VU 

where a = \/2(t V s) and «i = 

Let us now look at 1%. It is upperbounded by the factor i^o(2 3 / 4 ) 2 ai times 
the integral 

Since in the area of integration u\/ v = v = V ^ we are led to the product of 
2/5 by the one-dimensional integral 

a/2 (a - 2U) 5 / 2 

o 



U 
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Figure 5: The change of coordinates. 



where a 2 = dis- 
similarly I 3 is upperbounded by the product of A'Q(2 3 / 4 ) 3 aia2§§ times the 
integral 

a/2 ^- 2U ^ dUd V = a 3 a^ 



'U 

where a 3 = 12 3 ^ ■ One easily shows then that 
/„ < A^F(2 3 / 4 )"2 n (f[oi 



n"=i(2+3(i-i))y 

It can be verified by using a system for symbolic computation that < <x; < 1 
for all i > 1. One also notices that 

n on— 1 

H(2 + 3(j-l))>—(n-l)\, 

3=1 

therefore 



In < /^»(2 3 / 4 )"2"~ 1 3~ ( "~ 1) 1 

(n — 1). 



and this finishes the proof. 

□ 
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